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Abstract 

The  purpose  of  this  investigation  was  to  quantify  the 
effect  initial  structural  loading  of  a  graphite/epoxy 
member  has  upon  the  mass  removal  rate  during  laser 
ablation.  The  effective  heat  of  ablation  (Q*  =  enerqy 

absorbed/mass  removed)  was  used  as  a  measure  of  this 
efficiency.  A  simple  physical  model  of  the  important 
factors  affecting  the  graphite/epoxy  was  developed,  and 
predictions  were  made  of  the  effect  of  loading  on  Q* 

•  f  f 

A  three-dimensional  finite  difference  heat  transfer 
code  was  written  to  predict  the  temperature  distribution  in 
the  composite.  The  orthotropic  nature  of  the  thermal 
conductivity  tensor  in  the  plys  was  modeled,  to  accurately 
moael  the  heat  flow  from  the  irradiated  region. 

The  effect  of  thermal  and  mechanical  loads  upon  the 
stress  distribution  in  a  single  fiber  was  calculated,  and  a 
linear  decrease  in  Q*  „  with  increasing  stress  was 

•ff 

predicted.  The  coefficient  of  increasing  ablation 
efficiency,  (rj. ),  was  postulated  to  increase  linearly  with 
axial  stress  a  .  This  was  based  upon  the  hypothesis 
that  fracture  of  individual  fibers  will  be  a  process 
linearly  dependent  upon  applied  stress,  and  will  remove  a 
fraction  of  the  composite's  mass  without  requiring  the 
absorption  of  laser  energy. 
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Uniaxial  tensile  coupons  were  fabricated  from  Hercules 
AS4/3502  graphite/eDoxy  prepreg,  in  balanced,  symmetric 
laminates  laid  up  in  a  pattern  of  i 0/+60/-60 J  .  These 

n8 

specimens  were  place  under  tensile  loads  between  0%  and  50% 
of  the  laminate's  fracture  strength.  While  under  fixed 
grip  loading  conditions,  they  were  irradiated  with  a  10.6 
Ann  device  at  irradiances  between  5  and  26  kW/cm* . 

Linear  regression  of  data  taken  at  ~15kW/cmz  calculated 
17 .  =  0.013  kJ/gm  +  1 . 74  ( kJ/gnD/GPa**?  .  At  26  kW/cm*, 

a  x 

tk  =  0.087  kJ/am  +  3 . 36 ( k J/gm) /CPa*^  .  Therefore,  it  may 

d  x 

be  concluded  that  axial  fiber  stress  doe3  linearly 
decrease  the  10.6  Aim  laser  energy  needed  to  ablate  AS4/3502 
graphite/epoxy.  The  total  effect  seen  in  this  study 
was  less  than  20  percent  of  the  unloaded  value. 

It  is  recommended  that  further  analytic  modeling  of 
the  interaction  between  applied  load  and  fiber  fracture  be 
undertaken  to  allow  prediction  of  r>.  values  for  different 
materials  without  the  collection  of  experimental  data. 
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THE  EFFECT  OF  LOADING  ON  THE  LASER  ABLATION 
OF  GRAPHITE/EPOXY  COMPOSITE  MATERIAL 

Introduction 

The  laser-induced  ablation  of  fiber-reinforced 
composite  materials  has  been  studied  extensively. 

However/  the  majority  of  experiments  performed  have  been 
on  samples  of  material  without  any  externally  applied 
loads  (5/19).  Since  most  structures  of  interest  do  bear 
loads,  the  effect  of  these  loads  on  the  phenomenon  of 
ablation  is  important  in  understanding  the  interaction 
between  laser  radiation  and  structures. 

Although  load  enhancement  of  mass  removal  rates  has 
been  observed  experimentally  (17),  little,  if  any, 
analytical  modeling  of  laser  ablation  phenomena  has 
included  appropriate  mechanical  processes  to  explain  this 
observation  (19). 

Also,  global  structural  instability  has  typically  been 
modeled  as  simple  degradation  of  material  strength 
properties  at  elevated  temperatures  (19).  This  approach  is 
clearly  not  applicable  to  simple  mass  removal,  which  takes 
place  under  no  load  in  materials  needing  no  structural 
strength  ( 2 )  . 
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Examination  o£  solid  particles  ejected  from  ablating 
graphite/epoxy  composite  during  the  Very  High  Irradiance 
( HI - I )  and  Screening  Test  Series  (STS-1)  tests  (2,18) 
indicated  that  'rafts'  of  solid  laminae,  on  the  order  of 
1000-5000  microns  were  ejected  from  ablating  samples. 

This  was  under  conditions  of  no  pre-loading.  It  is 
conjectured  that  if  the  state  of  the  sblation  surface  is 
such  that  discrete  volumes  of  material  can  be  broken  off 
and  ejected  without  being  vaporized  by  laser  radiation, 
any  increased  stress  due  to  applied  loads  will  increase 
the  rate  of  material  fracture.  This  would  have  the  effect 

of  allowing  more  mass  to  be  removed  discretely,  and  lower 

* 

the  apparent  heat  of  ablation  (Q  -ff)- 
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II .  Theoretical  Analysis 


Temperature  Distribution 

The  prediction  o£  stress  states  in  the  sample  requires 
a  knowledge  of  the  temperature  distribution  created  in  the 
sample  due  to  laser  heating  as  well  as  the  loading 
distribution.  Various  one-  (ID)  and  two-dimensional  (2D) 
finite  difference  schemes  have  been  used  in  the  past 
(5,13-16).  However,  the  detailed  knowledge  of 
temperature  gradients  needed  to  predict  the  fracture 
mechanics  of  the  fibers  suggests  a  three-dimensional  (3D) 
model  is  needed.  The  anisotropic  nature  of  the  composite 
material  also  contributes  to  this  necessity. 

The  conservation  of  energy  on  a  infinitesimal  control 
volume  leads  to  the  differential  equation  of  heat 
conduction  in  a  solid  (10:5).  Left  in  terms  of  a  finite 
control  volume  it  is 

-  J  q  •  n  dA  +  J  g  dV  =  J  pCp  dV  (1) 

A  V  V 

where 

-  J  q  •  n  dA  «  heat  flux  through  boundaries  of  V  (2) 

A 

J  g  dV  =  heat  generation  in  V  (3) 
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J  pC  dT/*t  dV  =  rate  of  energy  storage 

p  within  volume  V  (4) 

The  heat  flux  term,  q,  for  an  isotropic  solid,  is 

proportional  to  the  temperature  gradient,  q  =  -kTT. 

For  the  anisotropic  case,  k  is  not  a  constant,  but  a  second 

order  tensor  (10:612) 


_  =  _  K 

ij  ix 


(5) 


A  discretization  of  the  volume  of  interest  into  small 
elements  (rectangular  prisms  were  used  for  simplicity  in 
this  case,  although  other  schemes  are  possible)  yields  an 
explicit  value  for  the  temperature  rise,  dT,  at  each 
element  during  a  small  interval  of  time,  dt.  This  rise  is 
added  to  the  existing  temperature  of  the  element,  and  then 
the  process  is  repeated  for  the  next  time  interval.  Phase 
changes,  ablation  and  temperature  dependent  properties  may 
be  incorporated  for  each  element. 

Neglecting  heat  generation  and  substituting  Equation 
(5)  into  Equation  (2),  and  approximating  derivatives  with 
first  order  terms,  yields  the  following  equations  for  the 
energy  flux  into  an  element  of  volume  V  . 

IwK 


*x 

l,J,K 
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V  Ax  J  1  1*^  X,J.K  X,J+i,  K J 


(6a) 
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where  the  zero  elements  o£  K  have  been  le£t  out. 

To  compute  the  temperatures  in  a  region  divided  into 
L  rows,  M  columns  and  N  layers,  an  (L+2)  x  (M+2)  x  (N+2) 
array  of  temperature  nodes  was  created.  Modes  on  the 
extremities  were  set  identically  equal  to  their  closest 
neighbors  on  the  inside.  This  ensured  zero  flux  (due  to 
zero  temperature  gradient)  on  all  boundaries,  creating 
insulated  boundary  conditions. 

The  only  energy  input  into  the  volume  came  from 
absorbed  laser  flux,  which  was  added  to  each  face  exposed 
to  the  irradiance. 

Like  ID  and  2D  explicit  schemes  the  size  of  time 
interval,  dt,  is  constrained  to  be  less  than  a  critical 
value  for  stability.  As  suggested  by  Torvik  (13:13),  a 
conservatively  stable  time  increment  was  chosen  by  using 
the  smallest  dt  necessary.  That  is,  the  smallest  distance 
increment  of  the  three  possible  was  used  to  calculate  dt. 

Existing  data  on  high  temperature  thermal  properties  of 
AS4/3502  are  extremely  scarce.  Data  for  similar  materials 
were  collected  and  applied  (3,5,8,19).  Thermal  property 
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data  on  AS4/3501-6  composite  is  much  more  prevalent,  and 
has  been  used  interchangeably  with  that  available  on 
AS4/3502  in  this  work.  The  computer  code  applies  these 
properties  as  a  function  of  temperature  through  the  point 
where  data  are  available.  For  temperatures  above  those 
where  data  exist,  the  last  recorded  value  of  the  property 
is  maintained. 

In  a  one  dimensional  fiber  reinforced  epoxy  laminae, 
the  principal  directions  of  thermal  conductivity  are 
obviously  along  the  fiber,  normal  to  the  fiber  in  the 
plane  of  the  laminae,  and  normal  to  the  plane  of  the 

laminae  (Figure  1).  The  conductivity  ,  ,  in  the 

* 

direction  along  the  fiber,  et ,  is  dominated  by  the 
conductivity  of  the  fiber  itself.  Data  for  this 
conductivity  are  approximated  by  the  High  Temperature 
Materials  Information  Analysis  Center  (HTMIAC)  (8).  The 
data  is  calculated  by  a  parallel  conductor  analysis  using 
13  W/m*°K  for  fiber  conductivity,  and  an  isotropic 
conductivity  value  for  the  matrix  resin. 

It  is  assumed  that  in  the  cured  laminate,  fiber  to 

A 

matrix  ratios  will  be  similar  in  the  laminae  plane,  e2,  and 

a 

normal  to  them,  e  .  Therefore,  laminate  conductivities  K 
and  Kaa  are  assumed  to  be  equal. 

In  the  layers  laid  at  ±  a  (60  degrees  in  this  case), 
the  principal  directions  no  longer  coincide  with  the  global 
axes,  and  the  conductivity  tensor  is  no  longer  diagonal. 


K  -  Kl  *60  1 


K  -  K[  -60 
K  •  K[  0°  ) 


mm 

jgssgsga 


k3,  e3 


k3i1 


K  =  K  a  II  I  oyer  s 


Figure  1.  Thermal  Conductivity  in  Orthotrcpic  Laminae 
The  conductivity  values  in  the  rotated  layers  are  found 
by  the  second  order  tensor  transformation  (11). 


K'  =  a.  .  a  .  K  , 

ij  v  k  j  l  kl 


where 


K'  is  the  tensor  in  the  rotated  axes  e',e',e' 

1 4  i  2 '  a 

/V  A  A 

K, ,  is  the  tensor  in  the  principal  axes  e  ,e  ,e 
kl  t'  z  a 

A  A 

a  . ,  is  the  direction  cosine  between  the  e.  and  e' 

V  k  k  i 


Letting  A.,  be  a  matrix  of  the  direction  cosines  a.., 

ik  ik' 

Equation  (7)  can  be  written: 


K'  =  A..K..A.. 

u  vk  kl  jl 


For  the  rotation  o£  the  sample  in  the  laminate 
described,  this  leads  to  two  non-diagonal  values: 


ii 


=  K. 


Si 


=  { K  -  K  )  cosa  sina 

it  21 


(9) 


The  only  conductivity  terms  calculated  in  the  computer 
code  used  £or  this  project  are  these  and  terms  on  the 
diagonal.  It  would  be  a  simple  matter  to  calculate  and 
include  the  remaining  terms  if  necessary,  to  provide  £or 
totally  anisotropic  constructions. 

The  material  properties  of  each  element  in  the  region 
are  calculated  £or  the  current  temperature.  The  flux 
romponents  are  then  estimated  from  the  current  temperature 
distribution,  using  Equations  (6).  Equation  (1)  is  then 
solved  for  the  incremental  temperature  change,  AT.  This 
Increase  in  temperature  is  then  added  to  the  cell's  current 
value,  and  the  region  is  checked  for  the  onset  of  phase 
change . 

This  phase  change  is  ablation  in  this  case,  but  it 
might  also  be  melting  or  sublimation  in  the  case  of  other 
materials.  When  the  phase  change  temperature  is  exceeded, 
the  amount  of  energy  absorbed  is  subtracted  from  the  cell's 
total  phase-change  energy,  and  the  temperature  is  set  back 
to  the  phase  change  temperature.  If  the  cell's  energy  of 
phase  change  is  exceeded  during  this  check  in  later 
iterations,  the  cell  is  'removed'  from  the  calculations  by 
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having  . ts  conductivities  multiplied  by  zero  in  the 
following  iterations. 

The  current  code  using  this  algorithm  is  listed  in 
Appendix  A.  It  remains  predominately  stable  by  the 
criterion  that  temperatures  monotonically  decrease  from 
the  irradiated  regions  toward  insulated  boundaries.  This 
is  the  case  for  all  but  one  or  two  cells,  which  remain  5 
or  6  degrees  below  the  initial  temperature  imposed  upon 
the  region;  but  no  adjacent  cells  suffer  this  anomaly,  and 
the  remainder  of  the  region's  temperature  values  increase 
as  expected  even  at  long  run  times  (high  iteration  count). 

No  further  investigation  was  made  into  this  difficulty 
since  the  code  agrees  within  5  percent  with  simple  ID  and 
closed-form  analytic  predictions  of  the  temperature  profile 
through  the  thickness  (Figure  2).  It  also  predicts 
experimental  results  (Figure  3)  and  the  reverse  temperature 
gradient  was  found  to  decrease  with  increasing  mesh  refinement. 

Stress  Distribution  in  a  Single  Fiber 

Once  the  temperature  distribution  has  been  calculated, 
the  mechanical  response  of  the  fibers  will  be  determined. 

It  is  expected  that  fracture  of  the  fibers  will  occur  at 
stresses  below  their  ultimate  strength  due  to  growth  of 
cracks.  Initial  cracks  in  the  fiber  exist  in  some 
distribution.  Thermal  expansion  stresses  may  cause 
fracture  of  some  fraction  of  the  fibers  by  forcing  unstable 
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(a)  Semi -Infinite  Solution  vs.  98  Node  ID  Code 


(b)  ID  Code  vs.  3D  Code 


Figure  2.  Comparison  of  Various  Finite  Difference 
Temperature  Predictions 
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growth  of  these  cracks.  This  is  assumed  because  of  the 
evidence  of  discrete  fiber  removal  from  various  composite 
materials  with  no  applied  loads  (2). 

To  calcu’ate  the  stress  growing  the  cracks,  o*t 
(Figure  4)  the  fibers  will  be  modeled  as  single  beams 
laying  on  an  elastic  foundation  (Figure  5).  This  is  to 
include  the  influence  of  the  resin  layer  beneath  the 
fibers . 

Modeling  of  Single  Graphite  Fiber 

Elementary  beam  theory  provides  the  following  equations 
for  stresses  in  a  prismatic  bar  (4:214): 

^ 

To  find  this  moment,  a  single  graphite  fiber  will  be 
modeled  as  a  round  Bernoull 1-Euler  beam  with  two  distinct 
regions.  One  region  extends  from  the  centerline  to  the  edge 
of  the  laser  beam.  The  next  region  extends  beyond  the  edge 
of  the  laser  beam  into  the  unheated  remainder.  This  region 
is  many  orders  of  magnitude  longer  than  the  fiber  diameter, 
d,  so  that  the  fiber  may  be  considered  to  extend 
semi-inf lnitely  (Figure  6). 

The  axial  stress  in  a  fiber  due  to  an  axial  load  is 
simply  the  load  divided  by  the  fiber  cross-sectional  area. 
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TEMPERATURE 
(degree  KELVIN) 


Elapsed  time  from  beam  on 
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Figure  3.  Comparison  of  MELT3D  Predictions  with 
Experimental  Data 
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Figure  4.  Single  Graphite  Fiber 
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Figure  5.  Beam  Element  on  Elastic  Foundation 


B 


w  -  Wg  V  —  Vg 
S  =  SB  M  =  M« 


Figure  6.  End  Conditions  for  Beams  AB  and  BC 


if  shearing  on  the  matrix/fiber  surface  interface  is 
ignored.  Vhile  this  shear  is  the  only  physical  mechanism 
for  transfering  load  to  the  fiber,  the  remainder  of  the 
analysis  will  be  concerned  with  the  heated  region  after  the 
matrix  may  be  considered  to  have  been  removed  by  pyrolysis. 
This  means  all  the  load  will  be  transferred  to  a  fiber  by 
interaction  with  the  matrix  in  the  region  away  from  the 
actual  laser  spot,  so  simple  uniaxial  stress  is  assumed. 

The  actual  value  of  this  stress  in  a  single  fiber  is 
found  by  calculating  the  total  strain  the  laminate  will 
exparience,  and  imposing  this  strain  on  a  single  fiber 

-x  -  W  *°  (11> 

where  c°  is  the  strain  of  the  mid-surface  (which  is  the 
strain  everywhere  in  the  cress  section,  because  balanced, 
symmetric  laminates  do  not  couple  axiai  strains  and 
midplane  bending). 

Although  the  stresses  in  the  direction  of  the  fiber, 
a ■  ,  will  actually  be  the  sum  of  axial  and  bending  moment 
components,  for  the  conditions  being  explored  in  this 
effort,  the  stresses  due  to  bending  are  Insignificant 
compared  to  stresses  due  to  axiai  loads.  See  Appendix  B 
for  details  of  this  analysis. 
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Stress  Distribution  in  the  Laminate  (9) 

In  a  given  ply  o£  the  laminate,  the  plane  stress 

(a  -  r  =  t  =  0 J  constltuit ive  law  can  be  written  as 
*  x*  y*  J 


where  the  reduced  stiffnesses  Q. .  are 


v  E 
u  i 


v 

IS  SI 


v  B 

SI  * 

1  -  V  V 

is  u 


and  for  a  lamina  at  an  angle  9  with  respect  to  the 

structural  axes,  the  stress-strain  relations  in  terms  of 

Ntransformed  reduced  stiffnesses"  Q. .  are 

w 
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in  which 


(13) 


Q  =  Q  cos4©  +  2  |Q  +  2Q  Isin*©  cos*©  +  Q  sin4© 

11  11  I  12  <5d  J  22 

Qia  =  |q44  +  Q2t  -  4Q^j  sin*©  cos*©  +  Q4j  £sin4©  +  cos4©J 
Q  =  Q  cos4©  +  2  [q  +  2Q  Isin*©  cos*©  +  Q  sin4© 

ZZ  It  I  AS  OO I  22 


Q,a  ®  COS’e  +(Q„-Q.,+2Q««)sln’e  cos  6 

=(Q1/01,-2Q„)5ln’9  cos  9  *  (o„-Q„*2a<Jsln  9  cos*9 

'■(Q,.+QM'2Q..-2Q«s)*ln‘®  cos's  +  ««.  I”1"*®  +  «»*«) 

The  strains  anywhere  through  the  thickness  of  the 
laminate  may  be  expressed  in  terms  of  the  mid-surface's 

strain  -j*0,  s° ,  y°  \  and  curvature  ,  at  ,  <t  l  such 
l  x'  y '  *y)  ^  *'  y  xyj 

that  for  any  given  lamina 


a 

x 

•  a 
y 

T 

xy 


Jk 


Z 


* 

< 

* 


at 

x 

at 

y 

at 

xy 


} 


J 
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Integrating  these  stresses  over  the  thickness  produces 
the  following  expressions  £or  force  per  unit  width  and 
moment  per  unit  width 


N  1 

r tyM 

a 

M 

f 

a 

X 

X 

r~ >  r  *. 

X 

N 

-  -  | 

O' 

.  dz  =  )  J 

a 

y 

J 

y 

Z-  J 

y 

N 

-t/* 

T 

k  =1  a 

T 

xy  J 

xy  J 

k-i 

xy 

where  zfc  and  are  de£ined  in  Figure  7.  This  can  be 

expressed  as 


where 


A.. 

*4 


I  ( O.  (■. 

k«t 
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In  (9)  it  is  shown  that  for  a  laminate  cured  at  a 
temperature  other  than  the  service  temperature,  'thermal 
forces  and  moments'  may  be  defined  by 


where  the  thermal  expansion  coefficients  for  an  arbitrary 
lamina  at  angle  0,  I®*#  ayr  a*y}  *  are  given  by  vector 

transformation  of  the  two  coefficents  in  the  lamina's 


coordinates. 


Finally,  by  defining  effective  forces  and  moments 
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oc 


ItHs-iHM 

we  can  invert  the  stiffness  matrices  A,B  and  D  and  express 
the  strains  and  curvatures  as 


With  the  layup  used,  effective  engineering  properties 
can  be  defined 


E  = 

X 


TTX 


ss 


where  h  is  total  laminate  thickness  such  that 


o 

c 

X 


O’  /B 


X  X 


(19) 


With  the  properties  of  the  laminae  used  here  being 
E±  =  20.8  x  104  psi  (1*3.0  GPa) 

E#  »  1.38  x  10*  psi  (9.51  GPa) 

^  -  0.31 

G  *  0.85  x  10*  psi  (5.86  GPa) 

=  1 . 42\ 
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it  follows  that  the  effective  modulus  of  the  total  laminate 
is  E  =  8.1  x  10*  psi  (55.8  GPa).  Therefore,  =  a  /E  is 

X  XXX 

the  mid-surface  strain  which  will  be  assumed  for  each 
fiber  in  the  laminate.  The  stress  in  any  lamina  will  be 
assumed  to  be  this  strain  multiplied  by  the  lamina's 
stiffness  in  the  x  direction.  This  is  a  simplifying 
assumption,  since  once  ablation  occurs  symmetry  and 
balance  are  destroyed,  creating  a  laminate  with  different 
stiffness  coefficients.  These  differences  will  be  assumed 
to  have  negligible  effect  upon  the  axial  stresses  in  the 
test  samples. 


Interaction  of  Thermal  and  Mechanical  Loads 


The  matrix  epoxy  absorbs  little  of  the  total  energy  of 
ablation,  and  practically  no  evidence  of  discrete  removal 
of  the  matrix  exists  (2).  Therefore,  the  discrete  mass 
removal  model  will  only  consider  fracture  and  discrete 
removal  of  the  fibers.  This  leads  to  a  formulation  of  the 
form 


« 


total 


V  * 

m  m 


*  vr*r 


(20) 


where  V  =  volume  fraction  of  matrix 
m 

V  3  volume  fraction  of  fibers 

But  *f  is  only  the  energy  absorbed  by  that  fraction  of  the 
fibers  which  ablate,  A  fraction  i?d  (3  1-t^  )  will 

be  removed  intact,  absorbing  essentially  no  energy. 
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Although  Q*ff  is  an  empirical  measure  of  all  the 
energy  delivered  to  a  material  divided  by  the  mass  removed 
during  irradiation,  it  is  a  reasonable  approximation  to 
the  actual  ablation  energy. 

From  Figure  A2,  Cp  is  less  than  2000  J/kg-K.  From 
300  °K  to  ablation  at  3800  °K,  the  energy  absorbed  will 
be : 

t  =  C  AT 

P 

=  2000  j^-^-  •  3500  °K 


■  7 


kJ 

gm 


*  20\  •  Q 


•rr 


The  total  energies  of  temperature  rise  to  the  ablation 
point  plus  the  ablation  phase  change  energy  will  both  be 
considered  in  Q*ff. 

So  we  have 

Projected  -  Reflected  -  Radiated  -  Conducted  * 

Energy  of  Temperature  Rise  ♦  Energy  of  Phase  Change 

Reflected  ♦  Radiated  +  Conducted  +  Energy  of 
m  Temperature  Rise  ♦  Energy  of  Phase  Change 

3  Material  Lost  as  Solid  +  Material  Lost  as  Vapor 

The  true  Q  would  be  the  ratio  of  the  last  two  terms. 

Looking  more  closely  at  Equation  (20),  the  energy 
absorbed  will  be  examined  on  a  layer  by  layer  basis.  Bach 
layer  is  assumed  to  be  either  entirely  fiber  or  entirely 
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matrix/  with  thicknesses  divided  in  the  same  ratio  as  the 
composite's  f iber-to-matrix  distribution.  The  energy 
absorbed  during  the  ablation  will  be 

*  .  =  7)  Q*Ajq  =  77  Q*PrK  Ah. 

aba  a  t  t  a  f  f  la**r  t 

where 

Q*  =  effective  enthalpy  of  ablation  of  fibers  alone 
Ahf  =  thickness  of  a  fiber  layer 
A.  =  A,  =  laser  spot  area 

laaar  l 

Ahf  -  thickness  of  a  fiber  layer 


and  the  other  symbols  have  their  previous  definitions.  In 
the  next  layer  of  matrix,  the  energy  absorbed  is  simply 


aba 


=  Q  Am  = 


Q*P  A.  Ah 

m  m  laitr  vn 


The  effective  Q*  calculated  during  the  ablation  of  all 
these  layers  would  be 


E  * 


<xbm 


E  An 


(21) 
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where  N  =  total  number  of  fiber  layers 
M  =  total  number  of  matrix  layers 

in 

NfAhf  =  hf/  total  fiber  thickness 
N  Ah  =  h  ,  total  matrix  thickness 

m  m  m 

and  note  that  the  total  thickness  ablated  is 

h  =  h  +  h 
t  m 

where 


h  =  h  V 

m  m 

hf  -  h  Vf  (22) 

Substituting  Equations  (22)  back  into  Equation  (21) 
results  in 


V,VrQf  +  (1-V*rVfQf  * 


‘•f  r 


P  v 

t 


P  v 

m  m 


*rvr  + 


(23) 


p 


where  p  is  total 
composite  density 


(i  -  »d)  *,o*  *  *.  q; 


(24) 


where  *f  ~  PfVf/p  38  i®88  fraction  of  fibers 
X  =  p  V  /p  =  mass  fraction  of  matrix 

m  m  m 
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becomes  simply 


It  can  be  seen  when  r?.  is  zero,  Q*  _ 

0  #Ii 

•  m  * 

Q  =  Q_*_  +  Q  x  I  which  is  the  limiting  case  that 

#  f  f  r  *  in  m 

would  have  been  arrived  at  by  simple  intuition. 

The  £orm  of  relationship  between  applied  stress  a  and 
discrete  removal  efficiency  7)d  is  suggested  by  the 
fracture  mechanics  of  a  single  fiber.  The  mass  of  discrete 
fibers  removed  will  be 


where 

N.  <=  number  of  discrete  fibers  removed 

d 

L  =  average  length  of  a  removed  fiber 

The  brittle  fracture  of  a  given  fiber  will  be 

governed  by  the  applied  stress  intensity  factor  (SIP)  Kx. 

When  this  applied  SIP  exceeds  the  critical  SIP  K  crack 

ic 

growth  will  occur.  If  additional  energy  is  available 

to  allow  K  to  remain  above  K  ,  unstable  crack  growth 
x  ic 

will  fracture  the  fiber  (1).  The  applied  SIP  is 
calculated  from 


K 


x 


(26) 


where  f  is  a  function  of  the  crack  geometry.  For  a  given 
constant  distribution  of  initial  cracks  of  length  a,  is 
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simply  proportional  to  applied  stress  .  if  damage  Is 
assumed  to  be  proportional  to  ,  this  suggests  a 
linear  variation  of  nd  with  the  applied  stress  a >  .  To 
account  for  the  possibility  of  some  threshold  stress,  a  , 
below  which  no  discrete  mass  removal  by  fracture  occurs,  a 
variation  will  be  assumed  of  -he  form 


d 


'dO 


a 

X 


(27) 


Available  data  on  this  phenomenon  is  too  sparse  to  be 
conclusive,  but  it  does  support  the  assumption  of  a 
monoton  leal ly  increasing  function  of  a  (17).  A  small 
amount  of  additional  data  collected  in  preparation  for  the 
primary  set  of  experiments  also  lends  credibility  to  the 
form  of  Equation  (27).  See  Appendix  D  for  details. 

Substituting  Equation  (27)  into  Equation  (24)  yields 

*  (i~',do  -  v.)  *f ♦  *„<c  <2» 

The  values  of  Qf  and  could  be  deduced  from 
experiments  conducted  at  o*  *  0,  or  assumed  to  be  (5) 

Q*  9g  43  kJ/gm 

Q*  9t  4  kJ/gm  (29) 

IW 

Using  the  mass  fractions  measured  for  the  particular 
specimens  used  in  this  thesis,  and  the  values  in  Equation 
(29),  the  predicted  effect  of  loading  on  Q*ff  will  have  a 
form  like  that  shown  in  Figure  8. 
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III.  Experimental  Apparatus 


Test  Specimens 

Fabrication  o£  the  test  specimens  is  detailed  in 
Appendix  C.  A  listing  o£  the  specimens  fabricated  and 
pertinent  pre-test  data  is  given  in  Table  I.  All 
specimens  were  laid  up  as  balanced,  symmetric  laminates. 

The  thinner  (0.244"  =  0.620  cm)  samples  were  (0,±  60) 

The  thicker  samples  (0.367"  *  0.932  cm)  were  (0,±  60) 

^  ixa 

Specimen  Instrumentation 

Type  K  ( chromel/alumel )  thermocouples  and/or  strain 
gauges  were  attached  to  some  specimens  (Figure  Clb).  All 
samples  were  weighed  be£ore  being  tested  to  determine 
their  mass. 

Tensile  Loading  Machine 

The  device  used  to  load  the  samples  in  tension  was  a 
hydraulically  actuated  £rame  built  by  the  Air  Force  Weapons 
Laboratory's  Directed  Energy  Weapon  E££ects  Branch 
(WL/TALE).  A  schematic  diagram  is  shown  in  Figure  9a,  and 
a  photograph  o£  the  grips  holding  a  sample  follows  in 
Figure  9b.  Loads  of  up  to  30000  pounds  of  force  (133400  N) 
were  applied  in  this  program.  The  force  was  measured 
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by  a  Strainsert  Universal  Flat  Load  Cell,  Model  FL50U-3DPKD, 
SN  06669-1  (50000  lb  capacity)  located  above  the  top  hyraulic 
grip  of  the  machine. 


TABLE  I 

Test  Specimen  Specifications 


Beam  Diameter  Spot  Average  Panel  Properties 


Shot! 

SamplelDff 

Minor 

(cm) 

Major 

(cm) 

Area 

(cm* ) 

Density  Resin 
(gm/cm*)  (Wt\) 

Voids 

(Wt%) 

1 

CH13989-1I5 

2.7 

2.9 

6.15 

NA 

NA 

NA 

2 

CH13989-1I3 

2.7 

2.9 

6.15 

NA 

NA 

NA 

3 

CH13989-4I1 

2.5 

2.7 

5.30 

NA 

NA 

NA 

4 

CH13989-4I2 

2.5 

2.8 

5.50 

NA 

NA 

NA 

5 

CH13989-4I3 

2.6 

2.8 

5.72 

NA 

NA 

NA 

6 

CH13989-4I5 

2.6 

2.8 

5.72 

NA 

NA 

NA 

7 

CH13989-5I1 

2.6 

2.8 

5.72 

NA 

NA 

NA 

8 

JH13189-1I3 

2.6 

2.8 

5.72 

1.59 

27.75 

1.32 

9 

JH13189-1I4 

2.6 

2.8 

5.72 

1.59 

27.75 

1.32 

10 

JH13189-2I3 

1.65 

1.8 

2.33 

1.59 

28.13 

0.92 

11 

JH13189-2I4 

1.65 

1.8 

2.33 

1.59 

28.13 

0.92 

12 

JH13189-2I5 

1.65 

1.8 

2.33 

1.59 

28.13 

0.92 

13 

JH13189-2I6 

1.65 

1.8 

2.33 

1.59 

28.13 

0.92 

14 

JH13189-2I7 

1.65 

1.8 

2.33 

1.59 

28.13 

0.92 

15 

JH13189-2I8 

1.65 

1.8 

2.33 

1.59 

28.13 

0.92 

16 

JH13189-1I3 

1.65 

1.8 

2.33 

1.59 

27.75 

1.32 

17 

CH13989-2I1 

1.2 

1.3 

1.23 

NA 

NA 

NA 

18 

CH13989-2I2 

1.2 

1.3 

1.23 

NA 

NA 

NA 

19 

JH13189-2I9 

1.2 

1.3 

1.23 

1.59 

28.13 

0.92 

20 

CH13989-3I1 

1.2 

1.3 

1.23 

NA 

NA 

NA 

21 

JH13189-1I4 

1.2 

1.3 

1.23 

1.59 

27.75 

1.32 

22 

CH13989-3»2 

1.2 

1.3 

1.23 

NA 

NA 

NA 

23 

CH13989-3I3 

1.2 

1.3 

1.23 

NA 

NA 

NA 

24 

CH13989-3I4 

1.2 

1.3 

1.23 

NA 

NA 

NA 

25 

CH13989-3V6 

1.2 

1.3 

1.23 

NA 

NA 

NA 

26 

JH13189-2I10 

1.2 

1.3 

1.23 

1.59 

28.13 

0.92 

27 

JH13189-1I6 

1.2 

1.3 

1.23 

1.59 

27.75 

1.32 

28 

CH13989-2#7 

1.5 

1.7 

2.00 

NA 

NA 

NA 

29 

JH13189-1I7 

1.5 

1.7 

2.00 

1.59 

27.75 

1.32 

30 

JH13189-1I5 

1.5 

1.7 

2.00 

1.59 

27.75 

1.32 

31 

CH13989-3I8 

1.5 

1.7 

2.00 

NA 

NA 

NA 

32 

CH13989-3I9 

1.5 

1.7 

2.00 

NA 

NA 

NA 

33 

CH13989-3I7 

1.5 

1.7 

2.00 

NA 

NA 

NA 

34 

CH13989-3I10 

1.5 

1.7 

2.00 

NA 

NA 

NA 

35 

JH13189-2I11 

1.5 

1.7 

2.00 

1.59 

28.13 

0.92 

36 

JH13189-2I12 

1.5 

1.7 

2.00 

1.59 

28.13 

0.92 
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Continuous  Wave  (CW)  10.6  urn  Laser  and  Associated  Diagnostics 

The  Electric  Discharge  Coaxial  Laser  II  (BDCL-II)  is  a 
10.6  Ann  infrared  carbon  dioxide  laser  capable  of  40  kV  of 
output  power.  As  shown  in  Figure  10,  in  the  test  program 
described  here  a  fraction  of  this  power  was  removed  by  a 
sodium  chloride  beam  splitter  to  be  measured  by  a  Coherent 
Model  213  power  meter  (SN  B0H375).  Diamond-lathed  copper 
mirrors  are  used  to  take  the  remaining  power  and  focus  it 
on  the  target  plane  in  the  shape  and  intensity  desired. 

At  the  target,  infrared  emission  from  the  target  is 
detected  by  a  Thermogage  Model  8000-1A  germanium  pyrometer 
(SN  3247).  This  data  provided  information  on  the  surface 
temperature  of  the  ablating  samples. 

16mm  motion  pictures  are  also  recorded  at  500  frames 
per  second  (fps),  using  a  Redlake  Hycam  II  Model  41-0064 
(SN333).  Video  tape  of  the  ablating  samples  was  recorded 
with  a  Sony  DXC-M3  Color  Video  Camera  and  a  Panasonic 
AG6200  VMS  tape  recorder.  Placement  of  the  camera  and 
recorder  made  recording  their  serial  numbers  impossible. 


j 


Shutter 

0 


Sp I i tter 


EDCL  II 


(  a  ) 


(b) 

Figure  10.  Laser  Optical  Train 
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IV.  Experimental  Procedure 


Specimen  Preparation 

The  sample  to  be  tested  Is  placed  In  the  upper 
hyraulic  grip  and  positioned  to  be  perpendicular  to  both 
the  axis  of  loading  and  the  Incident  laser  beam.  The 
lower  grip  Is  placed  around  the  bottom  of  the  sample  and 
any  Instrumentation  wires  are  routed  away  toward  their 
respective  recording  devices  (Pigure  9b).  1600  psi  (11 

MPa)  hyraulic  pressure  was  supplied  to  the  grips  to  secure 
the  specimen.  Data  recording  devices  are  checked  to  see 
that  proper  readings  are  being  received  from  each  of  the 
instruments . 

Laser  Optical  Train 

The  optical  train  is  adjusted  to  produce  a  desired 
spot  size  upon  the  target.  Accounting  for  the  expected 
device  output  power,  this  adjustment  controls  the 
irradlance  on  target.  Aberrations  in  the  irradiance 
profile  which  might  be  encountered  include  ellipticity  of 
an  otherwise  circular  beam,  irradlance  flucuations  at 
various  locations  within  the  spot  area,  and  power 
fluctuations  with  time  during  the  test. 
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The  £lrst  difficulty,  a  non-circular  beam,  occurs  when 
magnifying  mirrors  reflect  the  beam  at  a  finite  angle  with 
respect  to  their  centerline.  Since  it  is  Impossible  to 
reflect  the  laser  energy  perpendicular  to  the  mirror  and 
still  usefully  propagate  it,  some  ellipticity  must  remain 
in  the  beam,  unless  two  or  more  mirrors  can  be  used  to 
offset  each  others  effects.  Major  and  minor  diameters  are 
recorded  in  Table  I. 

The  next  problem  is  affected  by  the  ellipticity,  by 
alignment  of  the  laser  resonator  cavity,  and  by  Inherent 
design  of  the  EDCL-I1.  Both  ellipticity  and  spacial 
power  variations  can  be  quantified  by  measuring  the  burn 
profile  in  a  soft  material  like  plexiglass  before  testing 
the  sample,  or  to  a  less  accurate  degree,  by  the  shape  of 
the  sample's  burn  craters. 

The  last  problem,  temporal  fluctuations  in  power,  can 
only  be  measured  during  the  test  by  a  fast  response 
power  meter  measuring  scattered  radiation  from  some 
part  of  the  optical  train.  If  the  variations  seen  are 
significant,  then  they  must  be  accounted  for  after  the 
fact  in  the  analysis  of  the  test  data.  As  shown  in  a 
typical  plot  of  this  reading  (Figure  11),  there  were  no 
significant  fluctuations  during  these  tests. 
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Data  Collection  During  Irradiation 


Once  all  instruments  are  connected  and  operating 
properly,  recording  of  the  various  signals  is  begun,  the 
load  on  the  tensile  machine  13  raised  to  the  level  desired, 
and  the  laser  is  turned  on. 

After  1  to  2  seconds  are  allowed  for  the  laser  power 
output  to  stabilize,  opaque  shutters  keeping  the  beam  from 
going  to  the  target  are  removed.  The  time  allowed  for  the 
shutter  to  eclipse  the  beam  is  kept  well  below  5  percent 
of  the  total  irradiation  time  to  ensure  the  variation  in 
intensity  profile  is  insignificant  to  the  sample's  total 
interaction  with  the  radiation. 

After  a  predetermined  amount  of  time,  the  shutter  is 
returned  to  its  initial  position,  removing  the  radiation 
from  the  target.  The  laser  is  then  shut  down. 

Post-test  Data  Collection 

Any  significant  fibers  broken  off  during  the  test  which 
can  be  seen  downwind  of  the  purging  air  flow  are 
collected.  The  sample  is  then  weighed  to  determine  its 
mass  loss.  It  is  then  placed  in  a  plastic  bag  and  kept 
from  further  handling  until  photomicrographs  can  be  taken. 

After  the  photography  of  the  sample's  ablation 
surface,  the  depth  of  the  crater  ablated  in  the  sample  is 
measured  with  a  depth  gauge  or  other  micrometer. 
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V.  Experimental  Data 


Table  I!  and  Figures  12,  12,  ar.d  14  summarize  the 
results  obtained  on  the  EDCL-1I  device.  Sample 
calculations  showing  equations  used  to  reduce  the  data  are 
shorn  in  the  following  paraaraphs. 

Before  irradiation  tests  occurred,  two  samples  'A* 
and  ’ B *  were  pulled  in  tension  to  failure  in  the  tensile 
machine.  Sample  A's  load/strain  curve  is  shown  in  Figure 
12.  Sample  B  is  similar.  Compare  the  effective  modulus 
=  8.1  msi  (55.9  GPa)  to  that  obtained  on  similar  samples 
tested  on  an  AFIT  MTS  tensile  testing  machine  (also  ~8.1 
msi  (55.9  GPa)  from  Figure  C5).  This  is  in  excellent 
agreement  with  the  effective  modulus  calculated  in 
Chapter  II. 

The  strength  of  the  laminate  was  calculated  by  finding 
the  stress  which  places  the  ultimate  strain,  *  ,  on  the 

laminate 


«r  ,  .  •  E  «  =  8 .  Ixl0tfpsi*0 .0142 

ult  Hull 

=  115020  psi  =  793  MPa 
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0  .001  .002  .003  .004  .005  .006  .007  .008  .009  .01 

Strain  (in/in  :  mm/mm) 


■  Figure  12.  Stress  vs.  Strain  Curve  for  Test  Sample  A 

a 

The  load  which  applies  this  stress  to  the  laminate  is  the 
'ultimate  load'  used  to  normalize  the  data  in  the  Stress 
Ratio  column  of  Table  II.  For  0.247"  (0.627  cm)  thick 
samples,  it  is  approximately  284001bf  (126300  N).  For 
0.370"  (0.939  cm)  thick  samples,  425001bf  (189000  N). 

It  was  found  that  the  mass  losses  achieved  on  the 
samples  used  in  tests  1  to  9  were  too  small  to  allow 
accurate  measurements.  The  analytic  balances  available 
could  weigh  samples  less  than  160  grams  to  within  0.1 
milligram,  but  heavier  samples  had  to  be  weighed  on  a 

■  device  with  only  0.1  gram  resolution.  The  2"  (5.08  cm) 
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TABLE  IX. 

Laser  Interaction  Parameters  and  Results 


Shot 

• 

POWER 

on 

TARGET 

(kW) 

ZRRALIANCE 

fc) 

TENSILE 

LOAD 

UEE) — TlcFT 

STRESS 

RATIO 

LASER 

TIME 

(sec) 

MASS 

LOSS 

(gm) 

Q* 

mass 

m 

10 

31.77 

13.6 

0 

0 

0.00 

2.02 

1.95 

32.99 

11 

30.64 

13.1 

9400 

41.8 

0.22 

2.02 

1.98 

31.12 

12 

30.73 

13.2 

14800 

65.8 

0.35 

1.77 

1.82 

29.95 

13 

32.72 

14.0 

19000 

84.5 

0.44 

1.43 

1.56 

30.07 

14* 

31.68 

13.6 

18700 

83.1 

0.44 

1.61 

1.71 

29.90 

15 

31.60 

13.5 

14400 

64.1 

0.34 

1.76 

1.81 

30.65 

16 

31.60 

13.5 

9500 

42.3 

0.33 

1.51 

1.50 

31.89 

17 

31.77 

25.9 

0 

0 

0.00 

2.02 

2.28 

28.13 

18 

32.72 

26.7 

0 

0 

0.00 

1.27 

1.45 

28.64 

19 

31.25 

25.5 

5300 

23.6 

0.12 

1.02 

1.07 

29.84 

20 

31.51 

25.7 

4500 

20.0 

0.11 

1.01 

1.09 

29.30 

21 

31.68 

25.9 

5000 

22.2 

0.17 

0.63 

0.63 

31.83 

22 

31.68 

25.9 

9700 

43.2 

0.23 

1.03 

1.17 

27.99 

23 

31.51 

25.7 

9800 

43.6 

0.23 

1.11 

1.42 

24.61 

24 

31.25 

25.5 

14600 

64.9 

0.34 

1.24 

1.45 

26.72 

25* 

32.38 

26.4 

19000 

84.5 

0.44 

1.02 

1.45 

22. IS 

26 

32.03 

26.1 

15000 

66.7 

0.35 

1.02 

1.18 

27.59 

27 

31.60 

25.8 

4700 

20.9 

0.16 

0.78 

0.80 

30.96 

28 

31.68 

15.8 

0 

0 

0.00 

2.01 

1.95 

32.64 

29 

31.25 

15.6 

0 

0 

0.00 

1.27 

1.20 

33.15 

30 

31.51 

15.7 

2800 

12.5 

0.10 

1.01 

0.94 

33.68 

31 

30.99 

15.5 

2900 

12.9 

0.07 

1.76 

1.80 

30.16 

32 

31.25 

15.6 

5300 

23.6 

0.12 

1.51 

1.48 

31.92 

33 

32.29 

16.1 

7600 

33.8 

0.18 

1.52 

1.50 

32.74 

34 

31.77 

15.9 

9700 

43.2 

0.23 

1.28 

1.23 

32.93 

35 

31.60 

15.8 

11900 

52.9 

0.28 

1.01 

0.95 

33.56 

36 

31.77 

15.9 

14500 

64.5 

0.34 

1.01 

1.06 

30.09 

* 

Samples 

•14  and  125 

broke  in 

half 

dur ing 

testing 
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wide  samples  used  for  the  first  9  tests  were  all  greater 
than  160  grams.  The  0.1  gram  resolution  was  more  than  10% 
of  the  achieved  mass  losses.  This  poor  accuracy  makes  the 
data  unusable,  and  the  values  won't  be  reported  here. 

The  first  column  in  Table  II,  Power  on  Target  (kW)  was 
found  by  taking  power  out  of  the  EDCL-II  device,  measured 
as  described  in  Chapter  III,  and  plotted  on  a  laser 
printer.  The  laser  plot  was  then  redigitized  to  record  the 
final  value  representative  of  the  total  laser  power.  This 
number  was  multiplied  by  0.868.  This  number  was  suppled  by 
the  WL/TALE  personnel  as  the  factor  of  mirror  losses. 


P 


target 


P  .  * 

d»vie* 


0.868. 


This  number  has  an  uncertainty  of  ±  15%,  a  typical  value 
accepted  by  the  laser  effects  research  community  for 
ballistic  calorimeters.  The  next  column,  Irradiance 
(kW/cm*)  is  obtained  oy  measuring  the  major  and  minor  axes 
of  the  burn  spot,  with  an  uncertainty  of  ±  0.05  cm, 
calculating  the  burn  area,  and  dividing  this  area  into  the 
Power  on  Target. 


I 


torq«t 


n/ 4  (d  .  )  (d 

minor  mayor 


) 
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The  tensile  load  from  the  load  cell  is  also  plotted  on 
a  laser  printer  and  redigitized.  Definite  bounds  on  the 
uncertainty  associated  with  this  number  are  unavailable. 
The  stress  ratio  was  found  by  taking  the  load  and  dividing 
by  the  calculated  fracture  load  for  the  sample. 

Laser  time  on  target,  obtained  from  an  infrared 
detector  viewing  a  mirror  down  the  beam  train  from  the 
laser  shutter,  was  also  plotted  and  redigitized.  This 
uncertainty  is  no  larger  than  the  resolution  of  placing 
the  digitizer  cursor  on  the  plot.  The  uncertainty  in  time 
values  will  therefore  be  ignored. 

Mass  loss  values  came  from  weighing  the  samples  on  a 
analytic  balance  before  and  after  testing.  Uncertainty 
here  is  ±  0.2  milligrams.  This  is  less  than  0.06  percent 
for  all  samples  reported  and  will  also  be  ignored. 

The  effective  enthalpy  of  ablation  becomes 

P  •  t 

q*  s  Energy  =  target 

•ff  “  mass  loss  “  Am 

.  H-  *  I*  °»1 


Q*  =  15%  +  0%  +  0%  *  ±  15% 

•ff 

Although  this  uncertainty  is  larger  than  the  ordinate 
change  for  the  15  kW/cm*  data,  this  will  not  invalidate  the 
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conclusions  stated  in  Chapter  VI.  For  the  data  taken  at  26 
kW/cm2,  changes  larger  than  15  percent  are  seen. 

As  stated  in  Chapter  II,  the  effective  enthalpy  of 
ablation/  Q*  ,  is  based  upon  the  energy  delivered  to 
the  target.  Energy  lost  to  the  surroundings  will  not 
contribute  to  material  ablation,  making  the  value,  of 
Q*  higher  tha  ,  the  actual  enthalpy  value  of  the 
ablating  material.  The  most  significant  loss  mechanism 
for  material  ablating  at  3800°K  will  be  reradiation  of 
energy  away  from  the  target  and  conduction  away  from  the 
ablation  surface.  This  first  value  may  be  approximated  by 
calculating  the  energy  flux  of  a  body  at  temperature  T 
with  an  emissivity  c  (assumed  to  be  equal  to  the 
absorptivity  a) 

_,4  4 

q  =  c</Y  =  act T 


=  0.9 


5.67  x  10 


4 


where  a  is  the  Stefan-Boltzman  constant 
-  10.6  x  10**  V/m* 


=  1.06  kW/cm2  3s  8\ 


^13  kW/cm2} 


•»  4%  •  [26  kW/cm2} 
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The  second  loss,  conduction  away  from  the  ablation 
surface,  will  be  estimated  as  conduction  into  the  depth  of 
the  sample  and  radially  from  the  ply  being  ablated. 

Prom  Figure  B6,  the  average  thermal  gradient  through  the 
thickness  is  approximately  1.5  x  10s  °K/cm  at  the  ablation 
surface.  Prom  Figure  A2,  3s  1  W/m-K. 


‘cond 


=  K  7  T 


—  '  1.5  x  10*  ^ 
m°K  cm 


-1.5  kW/cm 


%  12%  •  ^13  kV/cm*j 

*6%  •  f*26  kW/cm*l 


Radial  conduction  will  also  remove  energy  from  the 
ablation  surface.  A  conservative  bound  on  this  loss  is 

estimated  by  taking  the  maximum  gradient  £l.5  x  109  °K/craJ, 

maximum  conductivity  (10  W/m-°K  from  Figure  A2)  and 
calculating  the  power  lost  in  the  ply  being  ablated.  For 
a  1cm  diameter  spot  and  0.0132  cm  thick  plys 

P  ^  ,  =  10  -^5-  •  1.5  x  10°  '  n  •  1  cm  •  0.0132  cm 

radial  m  *  K  CKI 

=  0.622  kW 
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Therefore,  the  reported  values  of  Q*#ff  are  too  high  by  a 
figure  on  the  order  of  22  percent,  compared  to  values 
computed  by  accounting  for  energy  losses. 

Figure  13  is  a  plot  of  the  data  for  26  kW/cm2. 

Figure  14  shows  13  kW/cm2  and  16  kW/cm2  data  on  the  same 
plot.  This  was  done  because  of  the  close  agreement  of  the 
unloaded  Q*  values  at  both  intensities.  The  magnitude  of 
the  correlation  coefficient  of  Figure  14,  being  higher 
than  that  of  the  26  kW/cm2  data  alone  (0.627  vs.  0.599), 
justifies  this  method  of  presentation. 

The  lines  fit  to  Figures  13  and  14  are  least  squares 
linear  regressions  of  Q*ff  on  stress  ratio  yielding  lines 
of  the  form 

<n  -  * ♦ »  ■  yr: 


Equations  (28)  and  (29)  are  solved  for  no  and  using  the 

values  for  A  and  B  above.  It  was  calculated  that 

a 


f)  +  V 
'o 


alt 


=  0.013  —  +  1.74  o' 

gro  <jfa 

for  data  taken  at  15kW/cm2. 

"a  -  °-087  |h  +  ^ •  36  !^|5  o- 

for  data  taken  at  26  kW/cm2. 

The  small  value  of  17  in  both  cases  raises  the 

O 

question  of  its  significance.  To  discern  this,  the 
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average  value  of  Q  taken  at  zero  load  was  subtracted 

•ff 

from  each  data  set.  A  straight  line  was  fit  through  the 
origin,  and  the  sums  of  the  residuals  squared  associated 
with  each  fit  were  compared  to  those  from  the  two 


parameter  slope,  B. 


For  15  kW/cm‘ 


kJ 


Q  =  32.96  ^  -  5.80 
•  ft  gm 


-  •  f— 1 

gm  KJ 


residuals  squared  =  18.2  a 
for  one  parameter  fit 

fa*  -  32.93  !£)  =-5.71  ^  •  f-^-1 

l  #ff  gmj  9®  1*^.1 

residuals  squared  =  23.0  a 
For  26  kW/cm2 

Q*  -  30.23  -  11.2  •  f-^-l 

•  ff  gm  gm 

residuals  squared  =  44.7 
and  for  one  parameter  fit 

K<  -  32  00  is)  ■  -17-30  is  • 


residuals  squared  =  55.7  a 

Fitting  a  one  parameter  fit  to  the  data  to  account 
for  an  offset,  r?o,  does  not  give  any  better  representation 
of  the  data  collected  here. 

Figures  15a  through  15g  show  representative 
photographs  of  the  test  samples.  In  addition  to  ply 
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separation  and  loose  fibers  hanging  from  the  ablation 
crater  surfaces,  global  structural  damage  is  obvious  at 
the  higher  load  levels.  This  type  of  damage  obviously 
contributed  (in  an  unknown  degree)  to  the  mass  loss 
enhancement.  Macroscopic  fiber  removal  occurred  during 
many  of  the  tests,  as  shown  by  fiber  bundles  seen  on  the 
floor  of  the  test  area.  A  typical  bundle  would  be  on  the 
order  of  2mm  wide  by  5mm  to  7mm  long.  No  attempt  was  made 
to  systematically  collect  these  bundles,  but  many  were 
found  downstream  of  the  nitrogen  sample  purge. 

Figures  16a  and  16b  show  scenes  from  the  video  tape 
taken  during  the  irradiation  of  sample  15.  The  darker 
bands  and  regions  are  areas  of  lower  temperature.  The 
orientation  of  these  bands,  which  coincide  with  ply 
orientation,  is  taken  as  evidence  of  discrete  ply  removal. 
A  solid  ply  breaking  off  and  being  carried  away  by  the 
purge  gas  would  expose  a  ply  beneath  it  at  lower 
temperature.  This  was  seen  in  MELT3D.F0R  predictions,  and 
has  been  observed  in  other  laser  effects  tests.  If  all 
plys  were  ablated  entirely,  without  discrete  removal,  new 
plys  will  not  be  exposed  until  the  layer  beneath  had 
heated  up  to  the  ablation  temperature,  and  no  darker  bands 
would  be  visible  on  film  or  video. 
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Figure  13.  Q*  versus  Stress  Ratio  for  26  kW/cm2 
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Figure 


14.  Q*  versus 


Stress  Ratio 


for  13  &  16  kW/cm 


(a)  Shot  110;  I  =  13  kW/cm*;  Stress  Ratio  =  0.0 
Figure  15.  Representative  Test  Sample  Photographs 
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(e)  Shot  113;  Close  Up 

Figure  15.  Representative  Test  Sample  Photographs  (Cont'd) 
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VI .  Conclusions  and  Recommendations 

As  Figures  13  and  14  show,  there  is  a  slight  trend 
toward  decreasing  Q*  at  higher  stresses.  The  data  indicate 
that  at  stress  ratios  explored,  only  a  9%  decrease  in  Q* 
was  seen  at  13-16  kW/cm*  and  19. 8\  decrease  at  26  kW/cm*. 

This  shows  that  the  initial  postulate  is  essentially 
correct.  That  is,  discrete  £Jber  removal  occurs  during 
ablation.  The  removal  by  a  fracture  mechanism  which 
appears  to  be  a  linear  function  of  applied  stress  (or  some 
linear  mechanism)  also  appears  to  be  verified  by  the 
experimental  data  collected. 

Therefore  it  is  concluded  that  graphite/epoxy  will 
ablate  at  rates  that  increase  linearly  with  applied  ixlal 
stress.  This  will  continue  until  the  applied  load  is  so 
high  that  the  tensile  specimen  will  fracture  due  to  the 
creation  of  an  initial  flaw  by  burning  off  one  or  two  plys 
on  the  surface  (12).  One  sample  did  this  during  the  tests 
at  LHMEL-1  (Appendix  D)  and  two  samples  broke  midway  during 
Irradiation  (likely  at  the  point  when  the  cross  section 
became  too  small  to  carry  the  applied  load)  in  the  EDCL-II 
test  series. 

No  specific  experimental  efforts  in  this  area  are 
recommended  for  the  immediate  future.  For  the  stress 
ratios  seen  here,  ablation  efficiency  increase  is  small 
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enough  that  it  may  neglected  for  many  engineering  purposes, 
and  the  unloaded  ablation  enthalpy  values  may  be  used. 

More  sophisticated  analytic  modeling  might  be 
undertaken  to  explore  what  factors  might  be  affecting  the 
particular  size  of  broken  fibers.  The  work  might  be  yield 
some  insight  to  the  rise  in  Q*  seer,  at  16  kW/cmz  between 

9  MI 

a/a  t  =  0.12  and  0.28. 

ult 

It  might  also  prove  useful  to  explore  the  change  in 
stiffness  of  a  laminate  as  a  cylindrical  hole  is  ablated 
from  the  sample.  Integrating  the  mass  loss  as  a  function 
of  instantaneous  stress  (instead  of  assuming  constant 
stress)  is  a  likely  next  step. 

If  any  further  experiments  are  conducted,  collection  of 
the  ablated  particles  must  be  attempted.  Differing  ply 
layups  may  also  be  Interesting  to  explore.  More  care  needs 
to  be  taken  when  designing  test  samples;  as  mentioned  in 
Chapter  IV,  analytic  balances  capable  of  weighing  larger 
masses  are  more  coarse.  It  is  imperative  to  build  samples 
light  enough  to  have  their  mass  loss  weighed  to  high 
precision.  A  conflicting  desire,  unfortunately,  is  the 
recommendat ion  to  test  samples  much  wider  than  the  laser 
beam  diameter,  to  eliminate  effects  of  the  laminate's 
edges . 

More  sophisticated  modeling  of  the  boundary  conditions 
on  a  fiber  might  be  undertaken;  specifically  the  effect,  if 
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any,  o£  layers  above  the  plane  o£  the  ablation  surface 

upon  the  fracture  characteristics  of  the  fiber.  Also,  a 

less  empirical  model  of  r>.  might  be  found  by  measuring 

<1 

the  distribution  of  initial  cracks  in  the  virgin  fiber,  and 
using  that  distribution  to  account  for  different  lengths  of 


removed  fibers. 


Appendix  A:  MELT3D.F0R,  Three  Dimensional  Heat 


Conduction  and  Ablation  Program 

This  appendix  shows  the  flow  of  logic  foe  this  program 
in  Figure  Al.  A  listing  of  the  code  follows.  Graphs 
showing  the  temperature  dependance  of  thermal  properties 
computed  in  this  program  are  given  in  Figure  A2.  A  listing 
of  the  subroutines  called  by  MELT3D . FOR  concludes  the 
appendix. 
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Check  interior  nodes 
for  pyrolysis  onset. 
Increment  ablation 
energy.  reset  temp 
t o  pyr o lysis  temp . 

IP  ob I ot 1  on  energy 
reached,  remove  cell. 


Pigure  Al.  Flowchart  for  MELT 3D. FOR 


58 


noon 


PROGRAM  MELT 3D 


C! 

C 

C 

C 

C 

C 

C 

c 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 


LAST  BDIT  17  October  1989,  21:40 
Created  2  August  1989  by  Capt  Joseph  L.  Hamrick, I I  GAE-89D 

This  is  a  3  dimensional  transient  heat  conduction  solving 
program  created  to  solve  for  the  temperature  distribution  in  an 
anisotropic  solid  by  finite  summation  of  the  heat  balance  on 
a  finite  control  volume.  The  differential 
equation  of  heat  conduction  is  left  in  integral  form  for  a 
differential  element  of  volume.  The  heat  fluxes  through  each 
face  of  the  cubic  elements  are  summed  algebraically,  and  the  net 
rate  of  increase  of  heat  in  the  element  (  in  Watts)  is  equated  to 
the  right  hand  side  of  the  equation,  (density*heat  capacity*rate 
of  temperature  change*volume  of  element;  or  rho  *  specheat  * 
dtemp/dt  4  vol).  This  is  solved  for  the  temperature  rise,  dtemp, 
and  added  to  the  existing  temperature  at  one  reference  node  of 
the  cube. 

The  heat  fluxes  are  discretized  in  terms  of  the  differences 
in  temperatures  at  discrete  nodal  points,  and  for  a  non-diagonal 
conductivity  tensor,  the  flux  through  any  side  is  simply  the  sum 
of  the  proper  tensor  terms  and  appropriate  temperature 
differences. 

This  program  is  only  intended  to  solve  for  the  instance  of 
prescribe  laser  flux  on  the  'top*  surface  (z*0)  and  INSULATED 
boundarys  everywhere  else.  This  is  justified  from  references  in 
the  thesis  which  show  convection  and  reradlation  to  be 
negllbible  for  the  interaction  regimes  of  interest  here. 


/  /  1 


qx(l,  j,k)»»«>  1  1  s*»=as>qx(  i+1,  j,k) 

1  qy(l,j+l,k)  1  / 

1  //  1  / 

1  //  1  / 

1  \//_  1  / 

1 _ 1/ 


THE  VARIABLES  USED  IN  THIS  PROGRAM  ARE: 


C 

C  LOOP  *  CRITBRIA  TO  RBAD  DATA  INTERACTIVELY  OR  FROM  A  FILE 

C 

C  TEMP ( I , J , K )  =  TEMPBRATURE  AT  NODE  (I,J,K);  APPROXIMATELY 

C  EQUAL  TO  AVERAGE  TBMP  IN  CELL  I,J,K  (deg  Kelvin) 

C  QX( I , J,K )  *  FLUX  IN  POSITIVE  X  DIRECTION  THROUGH  THE  FRONT  FACE 

C  OF  CUBE  I,J,K  (Watts) 
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C  QY( I , J,K )  *  FLUX  IN  POSITIVE  Y  DIRECTION  THRU  THE  LEFT  SIDE  FACU 
C  OF  CUBE  I,J,K  (VattS) 

C  QZ( I, J,K)  =  FLUX  IN  POSITIVE  Z  DIRECTION  (DOWN)  THRU  TOP  FACE  OF 
C  CUBE  I,  J,K  (Vatts) 

C  FLUX ( I , J , K )  =  LASER  FLUX  INTO  TOP  FACB  OF  CUBES  EXPOSED  TO  THB 
C  BEAN  <»at*-s;  (=  IRRAD(Ir J,K)*DX*DY] ) 

C  RK11( I, J,K, )  =  CONDUCTIVITY  TENSOR  ELEMENT  1,1  AT  I,J,K, (W/m*K) 

C  RK22(I, J,K, )  =  CONDUCTIVITY  TENSOR  ELEMBNT  2,2  AT  I,J,K, (V/m*K) 

C  RK12(I, J,K, )  =  CONDUCTIVITY  TENSOR  ELEMENT  1,2  AT  I, J,K, (W/m*K) 

C  RHO(I, J,K)  =  DENSITY  OF  BULK  MATBRIAL,  kg/m**3,  IN  CELL  I,J,K 

C  SPECHBAT ( I , J , K )  *  HEAT  CAPACITY,  W*s/kg*K,  IN  CELL  I,J,K 

C 

C  XMAX  =  LENGTH  OF  REGION,  CM 
C  YMAX  -  VIDTH  OF  RBGION,  cm 
C  ZMAX  -  THICKNESS  OF  RBGION,  cm 

C  NUMX  =  NUMBER  OF  CELLS  IN  X  DIRECTION 

C  NUMY  =  "  "  "  "  Y  • 

C  NUMZ  =  "  "  "  "  Z  DIRECTION,  =  (I  PLYS  ) 

C 

C  DX,DY,DZ  *  LENGTHS  OF  INCREMENTAL  DISTANCES,  meters 
C  VOL  =  DX*DY*DZ  ,  I.E.  VOLUME  OF  INCREMENTAL  CELL 

C 

C  COND?? (TBMP( I , J,K ) ,RK? )  =  SUBROUTINES  TO 

C  COMPUTE  TEMPERATURE  VARIATION  IM  PRINCIPLE  CONDUCUTIVITIBS 

C  RK1  AND  RK2 

C  DENS ( TEMP, RHO)  *  SUBROUTINE  TO  COMPUTE  RHO(TEMP) 

C  HEATCP( TEMP, SPECHBAT)  =  "  ■  "  SPECHBAT • TEMP ) 

C  MAKEFLUX  =  SUBROUTINE  TO  COMPUTE  IRRAD  VALUES  AND  EXPOSED  CELLS 

C 

C*ttt«tt*t**t*t*ttt«t**tt*ttt*******t*t*tttt***t*t**fttt***t***t*tt*t**t 

INTEGER  FLAG(10,10,10) 

CHARACTER*70  LABTXT 
CHARACTER *10  OUTFILE 
REAL  IRRAD 

C 

DIMENSION  TEMP (10, 10, 10),  FLUX(10,10,10),  RK11(10,10,10) 

DIMENSION  RK12(10,10,10),RK22(10,10,10),NSTACK(100) 

DIMENSION  RHO(10,10,10),  SPECKEAT(10,10,10) 

DIMENSION  SLUSH(10,10,10),  QX(10,10,10) 

DIMENSION  QY(10,10,10),QZ(10,10,10),RK33(10,10,10) 

C 

COMMON  /DATA1/  IRRAD, ABSORB,XMAX, YMAX, ZMAX, NUMX,NUMY, NUMZ, 

SRK 1MAX , RK 2MAX , RHOM I N , CPM I N , H ABL AT , T I MEMX , TI N I TL , S AFB 
COMMON  /DATA2/  RADIUS,  XFLUX,  YFLUX,  TPYRO  , ALPHA, ITYPE, IPRN, 

& I XOUT, I YOUT, I ZOUT, OUTFI LB, LABTXT 


BNTBR  DATA  DIALOG  TO  RBAD  LAST  RUNS  DATA  INPUT,  AND  ALLOW 
AMY  CHANGES  DESIRED 


€0 


ooooon  noonnnnonn  onnnnnonn  n  n  n  n  o 


C 


WRITB(*,4)  ' BMTBR  0  TO  CONTINUE  VITH  FILE  DATA,  1  TO  ENTER  CORREC 
AT  DATA  INTERACTIVELY* 

READlVdl)')  LOOP 


CALL  DATALK(LOOP) 
CALL  FLXTLK ( LOOP ) 


OPEN  FI LB  TO  OUTPUT  DATA  TO,  BLANK  FILE  NAME  FORCES  DOS  TO 
PROMPT  FOR  IT,  UNIX  NEEDS  BXPLICIT  NAME 


OPEN  (UNIT*12,FILB  *  OUTFILE  , STATUS* 'NEW ) 
WRITB(12, * ( A70) * )  LABTXT 


START  COMPUTING  INITIALLY  NBBDBD  QUANTITIES 
CHANGE  XHAX,  IN  CENTIMETERS,  TO  METERS 


DX  *  XMAX  /  ( 100 . 0 *FLOAT ( NUMX )  ) 
DY  =  YMAX  /  ( 100 . 0*FLOAT ( NUMY )  ) 
DZ  *  ZMAX  /  ( 100 . 0*FLOAT(NUMZ)  ) 
VOL  *  DX4DY4DZ 
PI  *  4. 000* AT AN (1.000) 

ALPHA  *  ALPHA4PI/180 .00 
NGONE  =  0 
T  =  0.00 


NOTE,  CASE  SPECIFIC  STATEMBMTS  FOLLOW: . 

THB  MAXIMUM  CONDUCTIVITY  FOR  STABILITY  PURPOSES  IS  BEING 
COMPUTED  FOR  THB  SPECIFIC  CASB  OF  K22  =  K33,  AND  ASSUMING  THB 
MAXIMUM  SUM  OF  K(ll)  IS  ALWAYS  LESS  THAN  THB  TRACB  OF  THB 
CONDUCTIVITY  MATRIX.  THEREFOR  ,  THE  MAXIMUM  IS  BEING  SET  TO 
THB  TRACE  OF  THB  MATRIX,  AMD  K22-K33. 


CONMAX  *  (RK1MAX  +  2.0*RK2MAX)/3.0 

DT  *  ( AMIN1 ( DX, DY, DZ ) **2 ) *RHOMIN*CPMIN/ ( 2 . 00*SAFBtCONMAX ) 
NTIME  -  IFIX  (TIMEMX/  DT  ) 


THB  FOLLOWING  SUBROUTINE  WILL  INITIALIZE  THE  STACKING  SBQUENCE  OF 
THE  PLYS,  TO  DICTATE  THE  ROTATION  OF  BACH  PLY  IN  THB  GLOBAL  COORD 
SYSTBM 


CALL  STACKR(NSTACK) 
C 
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C  INITIALIZE  THE  TEMPERATURE  AND  PROPERTIES  ARRAYS 
C  CHANGE  HABLAT,kJ/g»,  TO  J/kg  BY  MULTIPLYING  BY  10**6 

Ct*t*tt«t*****ttt**t*tt***t**tt**tt*t*t***t*****t********t*tt******t*t 

c 

DO  99,  I  =1,  NUMX+2 

DO  99,  J  =  1,  NUMY+2 
DO  99,  K  =  1,  NUMZ+2 

TEMP(I,J,K)  =  TINITL 
99  CONTINUE 

C 

DO  101,  I  =  2,  NUMX+1 
DO  102,  J  =  2,  NUMY+1 
DO  103,  K  =  2,  NUMZ+1 
C 

CALL  CONDI 1{  TEMP ( I , J , K ) ,  RK1) 

CALL  COND22(  TBMP(I,J,K),  RK2, I METAL) 

CALL  DENS  (  TEMP(I,J,K),  RHO(I,J,K)) 

CALL  HEATCP (  TEMP(I,J,K),  SPBCHEAT( I,J,K) ) 

SLUSH(I, J,K)  =  HABLAT*1 . 0E+06*RHO ( I , J , K ) *VOL 
FLAG ( I , J , K )  *  1 
C 

C************************************************************®******** 
C  NOV  ROTATE  THE  CONDUCTIVITY  TENSOR  TO  PLACE  BACH  LAYERS  PROPER 
C  TENSOR I AL  QUANTITIES  INTO  ITS  ARRAY  INITIALLY 

C 

CALL  ROTATE (  RK11(I,J,K),  RK22(I,J,K),  RK12(I,J,K),  RK33(I,J,K), 
t  NSTACX(K),  RK1,RK2,  ALPHA,  IMSTAL) 


103  CONTINUB 
102  CONTINUB 
101  CONTINUB 

C 

C** ******************************************** ****************** 

C****ttt******t**t«***tttt****t*tt**ttt*t***t**«**tt**t**t*t***** 

C  NBXT,  INITIALIZE  THE  FLUX  DISTRIBUTION  ON  THE  SAMPLE;  THIS 
C  IS  AG  AIM  CASE  SPECIFIC,  AND  IS  DETERMINED  BY  THE  SUBROUTINE 

C  LINKED  AS  'MAKBFLUX' ,  WHICH  GETS  ITS  VALUES  FROM  THE  COMMON 

C  BLOCKS 

c 

CALL  MAKBFLUX  (FLUX) 

C 

C**t**t*tttt**t**t**t*t**tt**tttttt*»ttttt»*ttt*ttt**ttt****tt**t 

C  NOV  ITS  TIME  TO  ACTUALLY  INCREMENT  THE  TIME  AMD  LOOP  THRU 

C  THB  REGION,  SUMMING  HEAT  BALANCES  ON  EACH  CUBB  DX*DY*DZ 

C 
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C  MOTE  1 1 ! ! 1 1 ! 

C  IP  AMT  INDEX  IS  1  or  ?MAX,  THE  TEMP  THERE  IS  SBT  EQUAL  TO 

C  INDEX-2  OR  ? MAX-1,  SO  THAT  NO  FLUX,  I.E.  INSULATED  BOUNDARY, 

C  IS  ENFORCED. 

C 

C  NOTE  ALSO,  THE  FLUXES  ARE  SET  TO  ZERO  IF  THE  NODE  CALCULATED 
C  HAS  ABLATED  OR  MELTED  AVAY,  OR  IF  THE  CELL  INTO  VHICH  IT'S 

C  HBAT  WOULD  PASS  IS  MISSING! 11!!! 

C***************************************************************** 


DO  200,  N  =  1,  NTIME 
T  =  FLOAT(N)  *  DT 

C************************** 

DO  201,  I  =  1,  NUMX+1 
DO  202,  J  =  1,  NUMY+1 
DO  203,  K  =  1,  NUMZ+1 
C 

QX( I, J,K )  =  FLOAT(FLAG( I , J,K ) *FLAG( 1+1, J,K) )*( (RK11( I , J, K ) * 
£  (TEMP( I , J , K )  -  TEMP(I+1,J,K))»DY*DZ/DX  )  +  (  RK12(I,J,K)* 

&  ( TEMP ( I , J , K )  -  TEMP( I, J+1,K) )  *  DZ  )) 

C 

QY( I, J,K)  =  FLOAT(FLAG(I, J,K)*FLAG(I, J+1,K) )*( (RK12 ( I , J,K ) * 
£  ( TEMP ( I ,  J , K )  -  TBMP( 1+1, J,K) )  ‘  DZ  )  +  (  RK22(I,J,K)  * 

£  ( TEMP ( I , J , K )  -  TEMP( I , J+1,K) )  *  DZ*DX/DY)) 

C 

QZ(I,J,K)  =  FLOAT(FLAG(I,J,K)*FLAG(I, J,K+1) )*(  RK33(I,J,K)» 
&  ( TEMP ( I , J , K )  -  TEMP(I,J,K+1))  *  DY*DX/DZ  ) 

C 

203  CONTINUE 
202  CONTINUE 
201  CONTINUE 
C 

C ********************************************************************* 

c ********************************************************************* 
c ********************************************************************* 
C  NOV  LOOP  OVER  THE  IMTERIOR  OF  THE  ENTIRE  MATERIAL  CUBE,  VHICH 
C  EXCLUDES  THE  PBRBPHBRIAL  NODES,  VHICH  ONLY  EXIST  TO  ENFORCE  THE 
C  BCs. 

C  THB  IMTERIOR,  HOVBVER,  SIMPLY  IS  THE  CONDUCTION  PROBLEM, 

C  UNLESS  CELLS  ARB  REMOVED  AND  THEN  LASER  FLUX  IS  APPLIED  TO  THB 
C  CELL  VHICH  USBD  TO  BB  BENBATH  IT. 

c****t******t**********************«*t****«*************************** 
C «••*•**•••*********•**•*•**•***•••***•**•**•**•*•****•*•*•*•****•**** 

c 

DO  210,  I  =  2,  NUMX+1 
DO  211,  J  *  2,  NUMY+1 
DO  212,  K  =  2,  NUMZ+1 
C 

BIN  *  FLUX ( I , J , K ) *DX*DY  +  QX(I-1,J,K)  -  QX(I,J,K)  +  QY(I,J-1,K) 

£  -QY ( I , J , K )  ♦  QZ( I, J,K-1)  -  QZ( I, J,K) 
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c 

DBLTAT  *  EIN*DT/  (RHOd,J,K)  *  SPECHEATd,  J,K)4  VOL) 

TEMPd,  J,K)  3  TEMPd, J,K)  +  DBLTAT 
C 

212  CONTINUB 
211  CONTINUB 
210  CONTINUB 
C 

C************************** ********************* ***•***************: 

c****************************************************************** 
C  NOV  GO  AROUND  EACH  BOUNDARIES  AND  SBT  THE  TBMPBRATURB  ON  BACH 
C  EDGE  NODB  EQUAL  TO  THE  TEMPERATURE  JUST  INSIDB  PROM  IT,  SO 

C  THAT  ON  SUCH  BOUNDARIES,  THERE  WILL  BE  NO  FLUX  THROUGH 

C  TO  THB  INTERIOR 

Ct  ********************************* ************* t*t**ttt*t********t: 

C 

DO  220,  I  =  1,  NUMX+2 
DO  221,  J  -  1,  NUMY+2 

TEMP  ( I ,  J ,  1 )  =  TEMP  ( I ,  J ,  2 ) 

TEMP( I, J, NUMZ+2 )  =  TEMP( I, J,NUMZ+1) 

221  CONTINUE 
220  CONTINUB 
C 

DO  230,  I  =  1,  NUMX+2 
DO  231,  K  3  1,  NUMZ+2 

TEMP ( I , 1 , K )  =  TEMPd, 2, K) 

TEMP(I,NUMY+2,K)  3  TEMP(I,NUMY+1,K) 

231  CONTINUE 
230  CONTINUE 
C 

DO  240,  J  ■  1,  NUMY+2 
DO  241,  K  3  1,  NUMZ+2 

TEMP(1, J,K)  3  TEMP(2, J,K) 

TEMP(NUMX+2, J,K )  3  TBMP(NUMX+1, J,K) 

241  CONTINUE 
240  CONTINUB 

C ****************************************************************** 

c *•**********•«***•***••**••*•*•***•*•**•**••*****•***••*•*«***•*•* 

c 

C  NOV  THAT  THB  NBV  TEMPERATURES  HAVE  BEEN  COMPUTBD,  RECALCULATE 
C  THB  PROPERTIES  IN  THB  REGION 

C  AGAIN,  IP  I METAL  IS  1,  SKIP  ANISOTROPIC  TENSOR,  SBT 
C  RK11*RK13RK223RK33  AND  RK12«0 

C* •***•*••***••***•••**«•*•«*••**••**•*•••*•**•****•*•**•••** ****** 
C* ft*************************************** ************************ 

DO  301,  I  3  2,  NUMX+1 
DO  302,  J  3  2,  NUMY+1 
DO  303,  K  3  2,  NUMZ+1 
C 

CALL  COND1K  TBMPd,  J,K) ,  RKl) 

CALL  COND22(  TEMP(I,J,K),  RK2, I METAL) 

CALL  DENS  (  TEMP(I,J,K),  RHO(I,J,K)) 
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CALL  HBATCP (  TEMP { I ,  J ,  K ) ,  SPECHBAT( I, J,K) ) 

C 

CALL  ROTATE!  RK11(I,J,K),  RK22{I,J,K),  RK12(I,J,K),  RK33(I,J,K), 

&  NSTACK(K ) ,  RK1,RK2,  ALPHA,  IMETAL) 

C 

C  IF  THIS  NODE  IS  NOT  AT  THE  PYROLYSIS  TEMPERATURE,  GO  ON  TO  THE 
C  NEXT  NODE.  IF  IT  IS  AT  PYROLYSIS  TEMPERATURB,  INCREMENT  THE 
C  'COUNTER'  KBEPING  TRACK  OF  THE  ABLATION  ENERGY  BEING  DUMPED  INTO 

C  THE  CELL,  AND  RESET  THB  TEMPERATURE  OF  THE  CBLL  BACK  DOWN  TO  THE 

C  PYROLYSIS  TEMPERATURE,  TPYRO. 

C 

IF  (TEMP(I,J,K)  .LE.  TPYRO)  GOTO  303 
C 

SLUSH(I, J,K)  =  SLUSH ( I, J,K)  -  RHO(I, J,K)*SPECHBAT(I, J,K)* 

& ( TEMP ( I , J , K ) -TPYRO ) * VOL 
TEMPI  I , J,K )  =  TPYRO 

C  IF  THE  'SLUSH*  ENERGY,  I.B.  ENERGY  TO  ABLATE  THE  WHOLE  CELL 

C  HASN'T  BEEN  REACHED,  GO  ON  TO  NEXT  CELL.  IF  IT  HAS,  REMOVE  IT 

C  PLACE  IT'S  FLUX  ON  THE  NBXT  CELL 

c. 

IF  (SLUSH(I,J,K)  .GT.  0.00)  GOTO  303 
FLAG ( I , J , K )  =  0 
FLUX! I , J,K+1 )  =  FLUX( I, J,K) 

FLUX ( I , J , K )  =  0.000 
TEMP ( I , J , K )  =  0.000 

Qtttttt **••** tt *************** **************** tttttttt ttttttttttttttttt 

C  INCREMENT  COUNTER  TO  SEE  HOW  MANY  CELLS  ARE  REMOVED 
NGONE  =  NGONE  +  1 

C  LINES  TO  OUTPUT  TO  FILE  AND  TO  THE  SCREBN 

C 

Qtttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttttt 

c 

WRITE! 12, 800 )  'CR',I,J,K,T 
WRITB(*,800)  'CR',I,J,K,T 
C 
C 

303  CONTINUE 
302  CONTINUE 
301  CONTINUE 

c 

IF  (  MOD  (N, IPRN)  .EQ.  0)  THEN 
WRITE!  12, '(  "TI '  ',F10.5) ' )  T 


65 


DO  310/  I  =  2,  NUMX+1,  IXOUT 

DO  310,  J  =  2,  NUMY+1,  IYOUT 

DO  310,  K  =  2,  NUMZ+1,  IZOUT 

C 

Y  =  10Q.0*DY*  FLOAT( J-2) 

X  *  100.0*DX»  FLOAT ( 1-2) 

Z  =  100.0*DZ*  FLOAT (K-2) 

WRITE(12, ' (1X/3F10.4/F8.0) ' )  X, Y,Z,TEMP( I,J,K) 

C 

310  CONTINUE 
END  IF 

c 

c 

IF  (  MOD  (N/ITYPE)  .EQ.  0)  THBN 
WRITB(*,803)  T 

WRITE(*,801)  ( 100 . 0*FLOAT( J-2 ) *DY,  J=NUMY+l,2,-IYOUT) 

DO  410,  K-2,  NUMZ+1, IZOUT 

VRITE( *, * )  'PLY  LAYER  I  '/K,'  DEPTH®  100.0*FLOAT(K-2)*DZ, 'CM' 
DO  411,  I  =  2,  NUMX,  IXOUT 
X  =  100.0*FLOAT(I-2)*DX 

WRITE(*/802)  X, (TBMP( I, J,K) ,  J  *  NUMY+1, 2, -IYOUT) 

411  CONTINUE 
410  CONTINUE 
END  IF 

Ctt**t*tt*t**tt**t***t*t***tt********t**ttt*t*t****t*ft**tt**tt*t£tt*t*« 

c 

c 

200  CONTINUE 
C 

CALL  DENS  (TINITL,  RHOMAX) 

DBLTAM  >  RHOMAX  *  VOL  •  NGONE  *  1.0E03 

WRITE( *, * )  'RHOMAX, VOL, NGONE  =  •  ,  RHOMAX, VOL, NGONE 

VRITB(*,*)  'HASS  LOSS  IS  ',  DELTAM 

C 

WRITE (12,*)  'NGONE  =  ', NGONE 
WRITBU2,*)  'MASS  REMOVED  ®  ',  DBLTAM 

800  FORMAT  (1X,A2,3I4,F10.5, 'CELL  IJK  REMOVED  AT  TIME  T') 

801  FORMAT ( IX, ' X  /  Y>',100F7.2) 

802  FORMAT(1X,F5.1,100F7.1) 

803  FORMATCl'/'TIME  IS  ',F10.6,'  SECONDS') 

STOP 

END 
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TEMP  (degress  KELVIN) 


Pigure  A2 .  Values  of  Thermal  Properties  Used  in  MELT3D.FOR 
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ooooon  n  n  n  n  n  o  n  n  n  n  n  n  n 


SUBROUTINE  CONDUCT,  RK1) 


SUBROUTINE  TO  COMPUTB  PRINCIPLE  (AXIAL)  CONDUCTIVITY,  Kll,  AS  A 
FUNCTION  OF  TEMPERATURE, T,  IN  DBGREES  KBLVIN,  BETWEEN  250X  AND 
3800K.  CONDUCTIVITY,  Kll  IS  RETURNED  IN  WATTS/M*K 


IF  (T  .GE.  2600)  THEN 
RK1  =  7.6 

ELSE  IF  (T  .GE.  1400)  THEN 

RK1  =  10.85  -  1.25E-03  *  T 

ELSE  IF  (T  .GE.  600)  THEN 

RK1  =  11.725  -  1.875E-03  *  T 

i3LSB  IF  (T  .GE.  250)  THEN 

RK1  =  -1.46101  +  T*(. 04390889-3. 982522E-05*T) 
ELSE  IF  (T  .LT.  250)  THEN 
RK1  =7.0 
I ERR  =  1 
END  IF 
RETURN 
BND 


SUBROUTINE  COND22  (T,  RK2,  METAL) 


SUBROUTINE  TO  COMPUTB  PRINCIPLE  (NORMAL)  CONDUCTIVITY,  K22,  AS  A 
FUNCTION  OF  TEMPBRATURE, T,  IN  DBGREES  KELVIN,  BETWEEN  250K  AND 
3800K.  CONDUCTIVITY,  K22  IS  RETURNED  IN  WATT3/M*K 


IF  (T  .GE.  2600)  THEN 
RK2  =0.4 

ELSE  IF  (T  .GE.  8001  THEN 

RK2  =  3.888888889E-02  ♦  1. 38888889B-04  *  T 
BLSE  IF  (T  .GB.  500)  THEN 

RK2  =  2.363333333  -  2.76666667E-03  *  T 

BLSE  IF  (T  .GB.  250)  THEN 
RK2  =  0.56  ♦  8.4B-04*T 
BLSE  IF  (T  .LT.  250)  THEN 
RK2  =  0.77 
I ERR  =  1 
END  IF 

METAL  *  0 

RETURN 

END 


SUBROUTINE  DBNS(T,P) 


SUBROUTINE  TO  COMPUTB  DENSITY,  P,  AS  A  FUNCTION  OF  TEMPERATURE, T, 
IN  DBGRBBS  KELVIN,  BETWEEN  250K  AND  3800K.  DENSITY,  P,  IS 
RETURNED  IN  kg/»**3 


IF  (T  .LB.  550)  THEN 
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P  =  1600.0 

ELSE  IF  (T  .GB.  800)  THEN 


P  =  1200.0 

ELSE  IF  (T  .GT.  550)  THEN 
P  =  2480.0  -  1.6  *  T 
END  IF 
RETURN 
END 

C*t*ttt*t*t*tt*t*t**t****t*ttt********tt*ttt**t****tttt«*tt***t**t*«*t* 

SUBROUTINE  HBATCP(T,CP) 

C*************************************************************ttt***ttt 

C  SUBROUTINE  TO  COMPUTE  HEAT  CAPACITY,  CP,  AS  A  FUNCTION  OF 

C  TEMPERATURE, T, IN  DEGREES  KELVIN,  BET VEEN  250K  AND  3800K. 

C  HEAT  CAPACITY,  CP  IS  RETURNED  IN  JOULES  /  kg  K 

IF  (T  .GB.  1500)  THEN 
CP  =  2000.0 

ELSE  IF  (T  .GB.  850)  THEN 

CP  *  741.825  +  T*(l. 484103  -  4. 316507E-04*T) 

ELSE  IF  (T  .GB.  500)  THEN 

CP  =  1131.4286  +  0.657142857  *  T 

ELSE  IF  (T  .LT.  500)  THEN 
CP  =  20.0  +  2.88  *  T 
END  IF 
RETURN 
END 

C*t***tt***t*****tt******t«t*tt*»ttt***t**t*tt*t**tt*ttt*t*ttt*t****«tt 

Ct**t**ttt*tttt**t*ttt***t*t*t**t*t**tt**t«t«***tttt*t*tt*t**tt**tt*tt* 

SUBROUTINE  STACKR ( NSTACK ) 

DIMENSION  NSTACK (100) 

COMMON  /DATA1/  I RR AD , ABSORB , XMAX , YMAX , ZMAX , NUMX , NUMY , NUMZ , 
&RK1MAX, RK2MAX,RHOMIN,CPMIN,HABLAT,TIMEMX, TIN ITL, SAFE 
C********************************************************************** 
C  FOR  A  (0,  +/-  ALPHA InS  LAYUP,  0  MEANS  0,  +1  MEANS  + ALPHA, 

C  -1  MEANS  -  ALPHA 

C********************************************************************** 
NPLYS  =  NUMZ 

DO  10,  I  =  2,  (NPLYS/2 )+l 

IF  (MOD (I -2, 3)  .EQ.  0)  THEN 
NSTACK ( I )  =  0 

ELSE  IF  ( MOD (1-2,3)  .EQ.  1  )  THBN 
NSTACK ( I )  =  1 

ELSB  IF  ( MOD (1—2,3)  .EQ.  2)  THEN 
NSTACK ( I )  =  -1 
END  IF 

10  CONTINUE 

C 

DO  20,  I  =  NPLYS/2  +  2,  NPLYS+1 

J  *  (I  -  (NPLYS/2  +1)  )*2  -  l 
NSTACK ( I )  *  NSTACK ( I -J) 

20  CONTINUE 

RETURN 
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END 

C*ttt*tt*t**t**t*ttt*tt«*t*tt**tttttt*****ttt**ttt*tttttt*ttt*t«tt«t**t 

SUBROUTINE  FLXTLK  (LOOP) 


CHARACTER* 70  LABTXT 
CHARACTER* 10  OUTIILE 

COMMON  /DATA1/  IRRAD, ABSORB, XMAX,YMAX,ZMAX,NUMX,NUMY,NUMZ, 
&RK1MAX, RK2MAX, RHOMIN,CPMIN, HABLAT,TIMEMX,TINITL, SAFE 

C 

COMMON  /DATA2/  RADIUS,  XFLUX,YFLUX,  TPYRO,  ALPHA, ITYPE, IPRN, 

4 I XOUT,  IYOUT,  IZOUT,OUTFILB, LABTXT 

C 

OPEN  (UNIT=11,FILE= ' 3DINPUT2.DAT', STATUS= 'OLD* ) 

READ  (11,*)  RADIUS,  XFLUX,  YFLUX,  TPYRO,  ALPHA, ITYPE, IPRN, I XOUT, 
&I YOUT, IZOUT 

READ  (11, • (A10,A70) • )  OUTFILE, LABTXT 
C 

DX  =  XMAX  /  ( 100 . 0*FLOAT ( NUMX )  ) 

DY  ••  YMAX  /  ( 100 . 0*FLOAT  ( NUMY )  ) 

DZ  =  ZMAX  /  (100.0 ‘FLOAT (NUMZ)  ) 

CONMAX  =  (RK1MAX  ♦  2.0*RK2MAX)/3.0 

DT  =  ( AMIN1( DX,DY,DZ ) *  *  2 ) *RHOMIN*CPMIN/ ( 2 . 00*SAFB* CONMAX) 

NT I ME  =  IFIX  (TIMEMX/  DT  ) 


IF  (LOOP  .BQ.  0)  THEN 


GO 

TO 

999 

ENb  xF 

CONTINUE 

WRITE( *, 

*) 

i  i 

WRITE( *, 

*) 

'THE  TIME  STEP  IS  ' ,DT, '  SECONDS' 

WRITEf  *, 

*) 

'BROKEN  INTO  ' ,NTIMB, '  INCREMENTS' 

WRITE**, 

*) 

•  1 

WRITE( *, 

*) 

1  f 

WRITE( *, 

*) 

'1 

RADIUS  OF  CIRCULAR  BEAM  IS  ', RADIUS,'  CENTIMETERS' 

WRITB(*, 

•) 

'2 

X  COORD  OF  BEAM  CENTBR  IS  ', XFLUX, '  CENTIMETERS' 

WRITE( *, 

*) 

'3 

Y  COORD  OF  BEAM  CENTER  IS  ', YFLUX, '  CENTIMBTERS' 

WRITB( *, 

*) 

'4 

ABLATION  OR  MELT  TEMP  IS  ', TPYRO, '  DEG  KELVIN' 

WRITE(*, 

*) 

'5 

♦  AND  -  WRAP  ANGLE  OF  LAYUP  IS  *, ALPHA, ’  DEGREES' 

WRIT8( *, 

*) 

'6 

INTEGBR  FREQUENCY  TO  TYPE  RESULTS  TO  CRT  IS',ITYPB 

WRITE( *, 

*) 

'7 

INTEGER  FREQUENCY  TO  SEND  RESULTS  TO  DISK  IS’, IPRN 

WRITB( *, 

*) 

'8 

DISK  FILE  TO  OUTPUT  TO  IS ( lowercase )  ',  OUTFILB 

WRITEf *, 

*) 

'9 

LABEL  TO  INSERT  AT  TOP  OF  DISKFILE  IS' 

WRITB( *, 

' ( A7  0 ) 

')  LABTXT 

WRITB( *, 

*) 

'10 

FREQUENCY  TO  SEND  X  DATA  TO  DISK',  I XOUT 

WRITB( *, 

*) 

'll 

FREQUENCY  TO  SEND  Y  DATA  TO  DISK',  IYOUT 

WRITBf  *, 

*) 

'12 

FREQUENCY  TO  SEND  Z  DATA  TO  DISK',  IZOUT 

WRITB( *, 

*) 

r  i 

WRITB( *, 

*) 

'Bnter  an  Integer  1  of  variable  to  change,  or  0  to 

icontlnue 

1 

WRITB( *, 

*) 

t  I 

WRIT8( *, 

*) 

1  t 
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READ!*,*)  NSVER 
C 

C  IF  NSVER  ISN'T  A  PERMISSIBLE  I  (1-9),  PROMPT  AGAIN 
C 

IF  ((NSVER  .LT.  1  .OR.  NSVER. GT. 12) .AND.  NSVsR.NE.O)  THEN 

VRITE( 4, 4 )  'REPLY  MUST  BE  0  OR  1  thru  12,  PLEASE  TRY  AGAIN' 

GOTO  017 
END  IF 

C 

GOTO  (1,2,3,4,5,6,7,8,9,10,11,12),  NSVER 

C* ****** 

C  IF  ANSVER  IS  0,  JUMP  OUT  OF  LOOP  AMD  VRITE  NEV  VARIABLES  TO  DATA 
C  FILE,  BY  JUMPING  TO  LINE  >100 

Qttttttttttttttttttttttt ********** tt**t**ttttttttisttttt**t***t*tttt*** 

GOTO  100 

001  VRITE(4,4) 'NEV  RADIUS,  CENTIMBTERS,  IS?' 

READ ( * , * )  RADIUS 
GOTO  017 

002  VRITE(4,4) 'NEV  X  COORD  OF  CENTBRLINE,  IN  CENTIMETERS,  IS?' 

READ( *, * )  XFLUX 
GOTO  017 

003  VRITE(4,4) 'NEV  Y  COORD  OF  CENTERLINE,  IN  CENTIMETERS,  IS?' 
READ(4,4)  YFLUX 
GOTO  017 

004  VRITE( *, * » *NEV  PYRO  OR  MELT  TEMP,  IN  KELVIN,  IS?' 

READ ( * , * )  TPYRO 
GOTO  017 

005  VRITE(4,4 ) 'NEV  VRAP  ANGLE,  DEGREES,  IS?' 

READ ( * , * )  ALPHA 
GOTO  017 

006  VRITE(4,4) 'INTEGER  FOR  SCREEN  OUTPUT  ?' 

READ ( * , * )  I TYPE 
GOTO  017 

007  VRITE(4,4) 'INTEGER  FOR  DI3KFILE  OUTPUT  ?' 

READ( 4, 4 )  IPRN 
GOTO  017 

008  VRITE( 4, * ) 1  NAME  FOR  DISK  FILE  ?' 

READ ( * , ' (A10) • )  OUTFILE 
GOTO  017 

009  VRITE( 4, 4 ) ’  LABEL  FOR  DESCRIBING  DISK  FILE  ?' 

READ(4, ' (A70) ' )  LABTXT 
GOTO  017 

010  VRITB( 4, 4 ) ' INTEGER  FOR  X  DISKFILB  OUTPUT  ?' 

READ(4,4)  IXOUT 
GOTO  017 

Oil  VRITB(4,4) 'INTEGER  FOR  Y  DISKFILE  OUTPUT  ?' 

RBAD(4,»)  IYOUT 
GOTO  017 

012  VRITE(4,4) 'INTEGER  FOR  Z  DISKFILE  OUTPUT  ?' 

READ!4,4)  IZOUT 
GOTO  017 
C 
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100  CONTI NUB 
REWIND  11 

WRITE! 11,*)  RADIUS,  XPLUX,  YFLUX,  TPYRO,  ALPHA,  I TYPE,  IPRN,IXOUT 
4, IYOUT, IZOUT 

WRITE! 11, ' (A10,A70) ' )  OUTFILE,LABTXT 
999  CONTINUE 
CLOSE  (11) 

RETURN 

END 

C***t*t«t*tt*****tt**t********tt**tt**ttt*t*t****ttttt****t*t**t******t 

SUBROUTINE  DATALK  (LOOP) 

REAL  IRRAD 

COMMON  /DATA1/  I RRAD , ABSORB , XMAX , YMAX , ZMAX , NUMX , NUMY , NUMZ , 
&RK1MAX,RK2MAX,RH0MIN,CPMIN, HABLAT, TIMEMX, TIN ITL, SAFE 

C 

OPEN  ( UNI T=10 , FILE= ' 3DINPUT . DAT  * , STATUS= ' OLD ' ) 

READ  (10,*)  IRRAD, ABSORB, XMAX, YMAX, ZMAX, HUMX, NUMY, NUMZ 

READ  (10,*)  RK1MAX,RK2MAX,RH0MIN,CPMIN,HABLAT,TIMEMX,TINITL,SAFB 

C 

IP  (LOOP  .EQ.  0)  THEN 
GO  TO  999 
END  IP 

C 

017  CONTINUE 

DX  =  XMAX  /  ( 100. 0*FLOAT( NUMX)  ) 

DY  =  YMAX  /  ( 100. 0*PLCAT( NUMY)  ) 

DZ  =  ZMAX  /  ( 100 . 0*FLOAT(NUMZ )  ) 

CONMAX  =  (RK1MAX  +  2.0*RK2MAX)  /  3.0 

DELTAT  =  ( AMIN1 ( DX, DY, DZ ) **2 ) *RHOMIN*CPMIN/ ( 2 . 0*SAFE*CONMAX) 


WRITE( *, * ) 
WRITE( *, * ) 
WRITE(*,*) 
WRITB( *, * ) 
WRITE!*,*) 
WRITE( *, * ) 
WRITE!*,*) 
WRITB( *, * ) 
WRITB( *, * ) 
WRITB( *, * ) 
WRITE( *, * ) 
WRITE!*,*) 
WRITB! *,*) 
WRITE!*,*) 
WRITE!*,*) 
WRITE!*,*) 
— 'ITS!  *,  *  ) 
WHITE!*,*) 
WRITB! *, * ) 
WRITE!*,*) 
WRITB!*,*) 


1  IRRADIANCE  INCIDENT  IS  ', IRRAD, '  W/c***2' 

2  SURFACE  ABSORBTIVITY  IS  ', ABSORB 

3  LENGTH ( X )  OP  SPECIMEN  IS  ',XMAX,'  CENTIMETERS' 

4  YMAX  (Y)  OP  SPECIMEN  IS  ' , YMAX,  '  CENTIMETERS' 

5  THICKNESS! Z)  SPECIMEN  IS  ',ZMAX,  '.CENTIMETERS' 

6  NUMBER  OP  X  NODES  IS  ',NUMX 

7  NUMBER  OP  Y  NODES  IS  ' ,NUMY 

8  NUMBER  OP  Z  NODES  IS  ',NUMZ 

9  MAX  AXIAL  CONDUCTIVITY  IS  ' ,RK1MAX, '  W/oK' 

10  MAX  NORMAL  CONDUCTIVITY  IS',RK2MAX,'  W/mK' 

11  MINIMUM  DENSITY  IS  ' ,RHOMIN, '  kg/n**3' 

12  MINIMUM  Cp  IS  ' ,CPMIN, '  J/  kg  K' 

13  THERMOCHEM  HEAT  OP  ABLATION  ', HABLAT, '  kJ/gm' 

14  MAXIMUM  TIME  DURATION  IS  ' ,TIMBHX, '  SECONDS' 

15  INITIAL  SAMPLE  TEMPERATURB  ' ,TINITL, ’  KELVIN' 

16  SAPETY  FACTOR  POR  TIME  INCREMENT' ,SAFB 
DT  will  t  ', DELTAT, '  Seconds' 

Enter  an  integer  I  of  variable  to  change,  or  0  to 


&continue ' 


WRITB!*,*)  '  ' 
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WRITE(*,*)  '  ' 

READ( *, * )  NSWER 

C 

C  IF  NSWER  ISN'T  A  PERMISSIBLE  »  (1-16),  PROMPT  AGAIN 
C 

IF  ((NSWER  .LT.  1  .OR.  NSWER. GT. 16) .AND.  NSWER. NE.O)  THEN 

WRITE{ *,*)  'REPLY  MUST  BE  0  OR  1  thru  16,  PLEASE  TRY  AGAIN' 
GOTO  017 
BND  IF 

GOTO  (1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16),  NSWER 

C* ************************************t******tt*************** ********* 

C  IF  ANSWER  IS  0,  JUMP  OUT  OF  LOOP  AND  WRITE  NEW  VARIABLES  TO  DATA 
C  FILE,  BY  JUMPING  TO  LINE  1100 

C************************************ ********************************** 

GOTO  100 

001  WRITB(*,*) 'NEW  IRRAD,  W/cm**2,  IS?' 

READ(*,*)  IRRAD 

GOTO  017 

002  WRITE( *, * ) 'NEW  SURFACB  ABSORBTIVITY  (0.0  to  1.0)?' 

RKAD( *, * )  ABSORB 
GOTO  017 

003  WRITE(*,*) 'NEW  LENGTH,  CM,  IS?' 

READ( *, * )  XMAX 
GOTO  017 

004  WRITB(*,*) 'NEW  WIDTH,  CM,  IS?' 

READ( *, * )  YMAX 
GOTO  017 

005  WRITB(*,») 'NBW  THICKNESS,  CM,  IS?' 

READ( *, * )  ZMAX 
GOTO  017 

006  WRITE(*,*) 'NBW  NUMBER  OF  X  NODES,  IS?' 

READ ( * , * )  NUMX 
GOTO  017 

007  WRITE(*,*) 'NEW  NUMBER  OF  Y  NODES,  IS?' 

READ( *, * )  NUMY 
GOTO  017 

008  WRITE( *, * ) 'NEW  NUMBER  OF  Z  NODBS,  IS?' 

READ ( * , * )  NUMZ 
GOTO  017 

009  WRITB(*,») 'NBW  MAXIMUM  AXIAL  CONDUCTIVITY,  W/  a  K,  IS?' 

READ ( 4 , * )  RK1MAX 
GOTO  017 

010  WRITB(*,*) 'NBW  MAXIMUM  NORMAL  CONDUCTIVITY,  W/  m  K,  IS?' 

RBAD(*,*)  RK2MAX 
GOTO  017 

011  WRITE(*,») 'NEW  MINIMUM  DENSITY,  kg/m**3,  IS?' 

READ(*, * )  RHOMIN 
GOTO  017 

012  WRITB(*,») 'NEW  MINIMUM  SPECIFIC  HBAT(Cp),  J/kg  K,  IS?' 

READ( *, * )  CPMIN 
GOTO  017 

013  WRITE (*,*) 'NEW  THERMOCHBM  HEAT  OF  ABLATION,  kJ  /  g«,  IS?' 
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READ ( * ,  * )  HABLAT 
GOTO  017 

014  WRITE(*,*) 'NEW  MAXIMUM  TIME, SECONDS,  IS?' 

READ(*, * )  TIMEMX 
GOTO  017 

015  WRITER,*) 'NEW  INITIAL  TEMPERATURE,  deg  KELVIN,  IS?' 

READ!*,*)  TINITL 
GOTO  017 

016  WRITE( *, * ) 'NEW  SAFETY  FACTOR  IS?' 

RBAD ( * , * )  SAFE 
GOTO  017 

100  CONTINUE 
REWIND  10 

WRITE( 10, * ) IRRAD, ABSORB, XHAX, YMAX, ZMAX, NUMX, NUMY, NUMZ 
WRITE(10,*)  RK1MAX,RK2MAX,RH0MIN,CPNIN, HABLAT, TIMEMX, TINITL, SAFE 
999  CONTINUE 
CLOSE  (10) 

RETURN 

END 

C**************************** ************************* t ***********! ***• 

SUBROUTINE  MAKEFLUX ( FLUX ) 

C****************************************************** **************** 

C  INITIALIZES  A  CIRCULAR  BBAM  OF  UNIFORM  INTENSITY,  IRRAD*ABSORB 
C  CENTERED  ON  THE  FACE  OF  RBGION 
£•***• *************** ***************************** ******** ************* 
CHARACTER* 70  LABTXT 
CHARACTER* 10  OUTFILE 
REAL  IRRAD 

DIMENSION  FLUX( 10, 10,10 ) 

COMMON  /DATA1/  IRRAD, ABSORB,XMAX, YMAX, ZMAX,NUMX,NUMY,NUMZ, 
&RK1MAX,RK2MAX,RH0MIN,CPMIN, HABLAT, TIMEMX, TINITL, SAFE 
COMMON  /DATA2/  RADIUS,  XFLUX,  YFLUX,  TPYRO,  ALPHA, I TYPE, I PRN, 

& I XOUT , I YOUT , I ZOUT , OUTF I LB , LABTXT 

C 

DO  100,  I  =  1,  NUMX+2 
DO  101,  J  =  1,  NUMY+2 
DO  102,  K  =  1,  NUMZ +2 
FLUX ( I , J , K )  =  0.0 
102  CONTINUE 

101  CONTINUE 
100  CONTI NUB 
C 

NXCL  =  NINT  (  FLOAT ( NUMX )  /  (XMAX/XFLUX) ) 

NYCL  =  NINT  (  FLOAT (NUMY)  /  ( YMAX/YFLUX) ) 

DX  =  XMAX  /  FLO AT ((NUMX)  ) 

DY  =  YMAX  /  FLOAT ((NUMY)  ) 

C 

DO  200,  I  =  1,  NUMX+2 
DO  201,  J  =  1,  NUMY+2 
C 

XSQUAR*(  DX  *  FLOAT (  IABS(NXCL-I ) )  )**2.0 
YSQUAR=(  DY  *  FLOAT (  IABS(NYCL-J ) )  )**2.0 
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PNTRAD=XSQUAR+YSQUAR 

DIFF=PNTRAD  -  (RADIUS*RADIUS  +  (0.5*(DX*DX  +  DY*DY)  )) 

C 

IF  (DIFF  .LT.  0.0)  THEN 

Ct«*tttt**t*ttt**«tt*4t**tt*ttt**tttt*tttt**t*ft*tt»t******fttt**tt****** 
C  CHANGE  FLUX  FROM  WATTS/CM* *2  TO  WATTS/M* *2  BY  MULTIPLYING  BY 
C  10000.0 

£t*t*«t*t***t****t*t***t**t*t*tt*tttt**t*tt****t*t***t*t**t*t****ttt*t* 

FLUX(T,J,2)=  IRRAD*ABSORB*10000 .0 
END  IF 
C 

201  CONTINUE 
200  CONTINUE 

C 

DO  300,  I  =  1,  HUMX 

WRITB( *, 900 )  (FLUX  (I,J,2),  J=1,NUMY) 

300  CONTINUE 
C 

900  FORMAT( IX, 100F4 . 0 ) 

RETURN 

END 

C**«*ttt*tt •»*****«***•*< •t***t***tt**t*t*ft*tt«tt«t*tt*tttt***t**«**t* 

SUBROUTINE  R0TATB(R11,R22,R12,R33,  NSTAK,  RK1,  RK2,  ALPHA, 
&IMBTAL) 

Q**t***t*tt*t«*t«tft*tttt*tt*tt***tt*t*t*t**t*t*ttt*t***ttt**t*t*ttt*tt* 

C  SUBROUTINE  TO  ROTATE  CONDUCTIVITIES  TO  HAKE  TENSOR  OUT  OF  THEM 
C  IF  IMETAL  IS  1,  SKIP  ANISOTROPIC  TENSOR,  SET  RK11=RK1=RK22=RK33 
C  AND  RK12=0 

C********tttt***t*t****t****t*t**t**t**t***t******t**t«t**tt**t*tt«t**t 

IF  (IMETAL  .EQ.  1}  THEN 
Rll  =  RK1 
R22  =  RK1 
R33  =  RK1 
R12  =  0.000 
GOTO  104 
END  IF 
C 

IF  (NSTAK  .BQ.  0)  THEN 
Rll  =  RK1 
R22  =  RK2 
R12  -  0.000 
R33  =  RK2 
C 

ELSE  IF  (NSTAK  .EQ.  1)  THEN 
R33  =  RK2 

Rll  =  RK 1 *COS ( ALPHA )*COS( ALPHA)  + 

&  RK2*SIN( ALPHA) *S IN (ALPHA) 

R22  =  RK2*COS( ALPHA) *COS (ALPHA)  + 

&  RK1*S IN ( ALPHA) *S IN ( ALPHA) 

R12  =  COS (ALPHA) *S IN (ALPHA) * (RK1-RK2 ) 

C 

ELSE  IF  (NSTAK  .EQ.  -1)  THEN 
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R33  =  RK2 

Rll  =  RK1*C0S (ALPHA) *COS (ALPHA)  + 

&  RK2*SIN( ALPHA) *3IN( ALPHA) 

R22  =  RK2*C0S( ALPHA) *COS( ALPHA)  + 

&  RK1*S IN ( ALPHA) *S IN (ALPHA) 

R12  =  COS (ALPHA) *SIN  ( -ALPHA) * ( RK1-RK2 ) 

END  IP 
C 

C* ************************************************** ******************* 

C  JUMP  DOWN  TO  HERE  IF  ISOTROPIC  MATERIAL 

C*****t***tt****«*t*t************************************************** 

104  CONTINUE 
RETURN 
END 
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Appendix  B:  Formulation  of  Beam  Equations  on  an  Elastic 


Foundation 

An  elastic  beam  in  flexure  on  a  foundation  with  a 
finite  modulus  is  modeled  as  a  simple  Bernoull i-Euler  beam 
element  subject  to  a  distributed  loading.  This  loading 
will  be  broken  into  two  distinct  parts:  p(x),  the  normal 
transverse  loading  component,  and  q(x),  the  reaction  of  the 
elastic  foundation.  If  the  foundation  is  assumed  to  be 
linearly  elastic,  q(x)  =  k-v(x).  The  differential  equation 
of  the  beam  in  Figure  B1  becomes 

El  w" "  +  kw  =  -  p ( x )  (Bl) 

The  homogenous  solution  to  this  equation  (7:4)  is 

w(  x)  =e^*(  A'  cos/?x  +  B' s  infix)  -  +  e-^*  ( C'  cos/?x  +  D' s  infix )  (B2) 
where  ft*  =  k/4EI . 

Without  any  loading  p(x),  the  portion  of  the  beam  from 
B  to  C  will  obey  the  homogenous  solution.  Since  the  fibers 
are  very  long  compared  to  to  their  diameter,  the  fiber  may 
be  considered  to  be  semi-infinite  in  length. 

The  requirement  for  the  displacement  and  slo4“>e  of  the 
fiber  to  vanish  at  large  x  (w  =  dw/dx  =  0;  x  ■>  oo) 
dictates  that  the  coefficients  of  must  vanish... 
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z,  w 


Figure  Bl.  Beam  Element  on  Elastic  Foundation 

wbc(x)  =  |<C'cosf?x  +  D'sin/?xJ  (B3a) 

w;c(x)  =  ne'f**  £(  -C'  +D'  )cos/?x  +  ( -C' -D' )sinf?xj  (B3b) 

w"(x)  =  2flze^*  f-D'cos^x  +  C'sin^x]  (B3c) 

®  c  v.  J 

W— (x)  =  2ftae~^x  [(C'+D')cos(3x  +  ( -C'  +D'  ) s inf?xj  (B3d) 

The  thermal  gradient  across  the  portion  o£  the  beam 
from  A  to  B  would  induce,  in  a  beam  element  entirely  free 
to  rotate,  a  curvature  of  1/Pt  where  is  the  radius  of 
curvature  (Figure  B2). 
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However,  the  section  AB  is  not  free  to  rotate,  but 
instead  has  a  curvature  (assumed  constant)  equal  to  the  sum 
of  thermal  curvature,  l/pt  plus  the  curvature  due  to  the 
bending  moment  Mb  applied  at  the  ends  of  section  AB 
(4:406).  If  the  matrix  resin  (elastic  foundation)  is 
removed  by  pyrolysis  at  700°K  and  the  fiber  remains 
intact,  at  temperatures  above  700°K  the  heated  sectir*.  AB 
may  be  considered  to  be  constrained  as  shown  in  Figu. e  B3, 
with  curvature  given  by 


A  B  C 


Figure  B3.  End  Conditions  for  Beams  AB  and  BC 
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W''(x)  =  ./*'  +  W'  *  . 

AB  homogtn«ou«  thermal 


Two  integrations  with  respect  to  x  yields 


WAB(X) 


*  f  JL  +  JL  ]  *  *1 

l  EI J  2 


+  Ax  +  B 


( B5 ) 


Symmetry  of  the  fiber  about  x-0  requires  that  w* (0)  =  0 
and  V(0)  =  Elw'''(0)  =  0.  Continuity  at  point  B  (x=r)  is 
required  of  the  deflection/  slope  and  moment.  The 
substitution  y  =  x-r  will  be  made  in  Equation  (B3)  so  that 
section  BC  may  be  thought  of  as  a  semi-infinite  beam 
starting  at  y=0.  This  yields  the  following  system  to  be 
solved . . . 


*  Y~  +  Ax  +  B 


( B6 ) 


w  (y)  =  e  ^  (C'cos/?y  +  D'sin/9y) 

B  G 


( B7  ) 


subject  to 


wi»(0)  =  0  (B8a) 

w' * '  (0)  =  0  (b) 


w  (x=r)  -  w  (y--0) 

AB  bc  * 


(c) 
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w'  ( x=r )  =  w'  (y=0) 

AB  BC ’ 


(d) 


M  -=  El  w"  (y=0) 

B  BC  J 


(e) 


where  M^,  A,  B,  C' ,  D'  are  to  be  determined  from  Equations  (B8) 


(  B8a  ) ,  (  B6  )  -*• 


( B8c ) 


(  B8d ) 


k  •  fe] 
(+  *  ")■ 51 


0  +  A  =  0-*A  =  0 


+  B  =  C' 


*  r  =  m-C'  ♦  D'  ) 


( B9  ) 


(BIO) 


-M_ 


(B8e)  ->  M  =  El  (-2^*D‘)  -♦  D'  *  - — 

2ftZEl 


(Bll) 


Since  w^'s  0,  Equation  (B8b)  is  solved  identically.  However, 
equlibrium  of  forces  in  the  2  direction  requires  that 
V(x=r)  *  0  =  V(y=0).  This  leads  to 


0  =  2ft  ( C*  +  D'  )  -►  C'  =  -D' 


( B12 ) 


(BIO), (Bll), (B12) 


(  k  *  E']r  ‘ 


2ftD' 


=  2ft 


fJL-1 

l  2ft*EI  J 
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( B13 ) 


r_  _  V 

Pt  "  ft El  El 


(  Bll ) ,  ( B13  )  -#■ 


rEI _ 

♦  r]2^*EI 


*  ') 

(  B12  ) ,  (  B14  )  -»•  C'  =»  -D* 


♦  r) 

( B9  ) ,  (  B13 ) ,  ( B15 )  -► 


( B14 ) 


( B15 ) 
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B 


V 


2 


221 


'  A  l?1  * r  ■  ik\ 


( B16  ) 


where  y  s  i  +  r 


Note  two  interesting  limits  on  these  constants; 


lim  M  = 
ft+ao  " 


-El 


P 


t 


which  is  the  solution  for  a  double-fixed  beam  subjected  to 
a  temperature  gradient  (4:406).  This  agrees  with  the 
intuitive  observation  that  a  beam  section  BC  on  an 
infinitely  stiff  foundation  will  act  like  a  fixed  wall  to 
the  section  AB. 


Also,  lim  y-z ,  which  implies  lim 

ft+aa 


B  = 


r 


(r  -  r  - 


This  also  agrees  with  the  solution  for  a  double  fixed  beam. 
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The  reduction  in  bending  moment  in  the  section  AB  can 
be  seen  most  easily  by  defining  a  quantity 


M*  s 


-Ei/pt  r 


(B17) 


N  =1  implies  the  bending  moment  in  AB  is  identical  to  that 
of  a  fixed  wall  beam.  The  effect  of  laser  beam  radius  r 
and  foundation  stiffness  ft  is  shown  in  Figure  B4,  which 
clearly  shows  that  for  spot  sizes  and  foundation 
stiffnesses  seen  in  the  experimental  setup,  M%1 . 


M*(j8,r) 


b'  / 


-  r=0. 10 

-  -  r=0.25 

—  r=0.50 

—  r=  1 .00 

I - 1 - 1  I  I - 1 - 1 - 1 

25  50  75  100  125  150  175  200 

fi  (cm-’) 


Figure  B4.  M  versus  ft  and  r 
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To  estimate  ft,  first  estimate  the  foundation  modulus 

k,  which  is  the  ratio  of  beam  loading  per  unit  length  to  the 

deflection  (7).  As  shown  in  Figure  B5,  for  a  length  of 

beam  L  with  force  p  on  it,  the  pressure  on  the  foundation 

of  modulus  E  will  be  a  =  — — .  This  will  also  be  the 
m  dL 

stress  in  the  foundation  which  causes  strain  c  =  6/ h, 
where  h  is  the  thickness  of  the  foundation. 

cr  6  p 

£  ~  E  “  h  =  d  L  E 

m  m 

Since  k  '  L2?'  we  get 


d  L  E 

k  =  E.  _ " 

L  ph 


d  E 


Assuming  a  foundation  of  epoxy  resin  as  thick  as  the 
distance  between  fibers  (assumed  equal  to  a  fiber 
diameter ) 

d  E 


k  =s 


'm  8xl04cm  •  2.4xl07  N  cm-2 


8x10  cm 


Therefore, 

<t  -  " 


4  E  I  4B  d4 


_ 2.4x107N  cm~2 _ 

4x2 . 35xl07Ncm~2  — ^8xl0-4cmj 


ft  as  1890  cm  >>200  cm" 


8b 


0 


.2 


1 


1.2  1.4 


~i - 1 - r 

.4  .6  .8 

Time  (seconds) 


Figure  B6 .  Temperature  Gradient  in  Graphite/Epoxy 
at  lOkV/cro2  via  MELT3D.FOR 
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( 


M_  % 


For  a  laser  beam  radius  of  r=1.0  cm  and  ft>>200  cm 


El/pt. 


The  maximum  axial  stress  from  M  will  be 


_  M  y  _  M(d/2 ) 
"  ‘  '  I 


X  I 

M 


.*.  a 


E  I  (d/2) 


x  p  I 

M  t 


( B18 ) 


From  MELT 3D. FOR,  the  magnitude  of  the  thermal  gradient  at 
10kW/cm2  is  <  200/000  °K/cm  (Figure  B6 ) .  Across  the  first 
surface  fiber ,  the  temperature  difference  is  AT  %  150°K. 


From  Chapter  III, 


f1  I  °  Tul 


Pt  a  AT 


( B19 ) 


8x10 


4cm|l- 


3.3xl0'7  °K~ 4  *  380 0°K 


-3.3x10  7  '  150°K 


=  -16.1  cm 

Therefore  from  Equation  (B18)  and  Equation  (B19) 


~  E  (d/2) 
~  P. 


-*  _ _ -*  8x10  cm 

2.4x10  Ncm  •  - 2 - 


-16. i  cm 


-2 


%  -600  N  cm  =  -6  MPa 
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I 


But  a  ,  =  3795  MPa,  so 

ult 


3795 


=  0.16%. 


But  a 


-  E  c  .  When  loaded  to  2/3  fracture  load,  c  =  — x- 
X  ,  1  t  '  4  3 

axial 


2  „ 

'»  ■  3  E. 

axial 


X  =  2  (-°.:0148)  235x10°  Pa 

"axial  3 


=2.31  GPa 


axial  2.31  GPa 


=  61% 


6MPa 

2310  MPa 


=  0.26% 


Therefore,  bending  stresses  are  appropriately  negligible 
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Appendix  C:  Fabrication  of  AS4/3502  Graphite/Epoxy 


Test  Specimens 


The  graphite/epoxy  (GE)  specimens  tested  in  this 
program  were  laid  up,  cured  and  cut  to  size  by  the 
Structures  Survivability/Supportability  Group  of  the  Flight 
Dynamics  Laboratory,  Wright  Research  and  Development  Center 
(WRDC/FIBCA) .  Specimen  specifications  described  in  Figure 
Cl  were  given  to  the  thesis  sponsor,  and  fabricated  as 
described  in  this  appendix. 

Layup  of  the  panels  was  performed  by  a  vacuum  bag 
process  as  shown  in  Figure  C2. 

After  the  panels  are  laid  up  and  sealed  under  vacuum 
(Fig  C2a)  they  are  cured  in  an  autoclave  in  a  manner 
similar  to  that  suggested  in  Reference  6:20  (Figure  C3). 

1.  Pull  20  in.  (508mm)  minimum  vacuum  on  part. 

2.  Place  in  autoclave. 

3.  Raise  temperature  to  275°F(135°C)  at  3- 
5°F/minute  ( 2-3°C/minute ) . 

4.  Hold  for  15  minutes  at  275°F(135°C)  under 
vacuum  pressure  only. 

5.  Pressurize  autoclave  to  85psi(586  kPa). 

6.  Hold  at  275°F(135°C),  85psi(586  kPa),  and 
vacuum  on  the  part  for  45  minutes. 

7.  Raise  temperature  to  350°F(177°C)  at  3- 
5°F/minute  ( 2-3°C/minute ) . 

8.  Hold  for  2  hours. 

9.  Cool  part  to  150°F(66°C)  in  not  less  than  45 
minutes,  maintain  pressure  and  vacuum. 

10.  Remove  from  autoclave  and  unbag. 
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TYPICAL  ACCEPTANCE  PANEL 
FABRICATION  SEQUENCE 


I.  Baseplate:  Aluminum  1/4 to  1/2  in.  thick 

2  Cork  dam:  Cork  1/8  x  1  in.  with  pressure-sensitive  adhesive  backing  (Corprene)  or  equivalent 

3  Release  film:  Teflon,  nonperforated  0.001-0  004  in.  thick 

4  Release  fabric.  Fabric  enfab  TX  10-40  release  (porous)  or  equivalent 

5.  Prepreg  layup 

6  Release  fabric:  Fabric  enfab  TX  10-40  release  (porous) 

7.  Resin  bleeders:  Cloth,  fiberglass  No  120  (as  per  calculation  on  following  page) 

8  Release  film:  Teflon,  nonperforated  0  001-01  004  in.  thick 

9  Caul  plate:  Aluminum,  0  030  in.  minimum 

10.  Tape:  Pressure-sensitive,  green  polyester  silicone 

II.  Air  bleeder  Style  1 581  glass  or  equivalent 

12  Vacuum  bag:  Film  capron  80,  hi-temp  nylon,  0  002  in  thick 
13.  High  temperature  sealant:  Schnee  Morehead 

Mold  release.  Frekote-33  or  equivalent 


(a)  Typical  Vacuum  Bag  Layup  (6:23) 
Figure  C2.  Fabrication  of  Test  Samples 
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(b)  Layup  Process 


(c)  Layup  Process  (Continued) 
Figure  C2 .  (Continued) 
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(d)  Layup  Process  (Concluded) 
Figure  C2 .  (Concluded) 


350  °F 


The  cured  panels  were  examined  by  computerized  axial 
tomography  (CAT)  scan  to  reveal  any  areas  of  incomplete 
interply  bonding.  Panels  fabricated  with  thermocouple  wire 
embedded  in  them  were  x-rayed  to  discern  any  shifting  of 
the  thermocouple  positions  that  may  have  occurred.  No 
flaws  were  visible  to  either  of  these  techniques. 

Evidence  that  the  structural  integrity  of  the  specimens 
was  as  high  as  practically  achievable  was  also  found  in 
post  mortem  inspection  of  e  failed  specimen.  One  sample 
was  tested  in  tension  to  failure  (Figure  C4).  At  87%  of 
the  predicted  strength,  failure  occured  at  the  grip  fixture 
(note  step  in  output  created  by  delay  in  switching  the 
amplifier  gain  as  strain  gauge  output  changed  scales). 


0  .002  .004  .006  .008  .010  .012  .014 

Strain  (in/ln  :  mm/mm) 


Figure  C4 .  Load  versus  Strain  Gage  Output  for 
Sample  JH13189-4  13 
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Failure  was  limited  to  the  immediate  area  of  the  fracture. 


Virtually  no  delamination  was  observed  except  in  the 
midplane  of  the  specimen,  where  the  failure  remained 
between  the  intially  separated  ply.  The  strong  bonding  of 
fibers  to  matrix  which  prevented  this  delamination  from 
running  into  other  plys,  and  the  nearly  Intact  fracture 
surface  show  that  the  lamination  process  in  this  specimen 
was  very  sucessful  (Figures  C5-C8). 

Lastly,  fiberglass  end  tabs  were  attached  to  specimens 
intended  for  tensile  testing,  and  panels  were  cut  to  size 
using  a  diamond  saw. 
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Sample  JH13189-4 
Surface  (20X) 
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#3, 


Delamination 


I 

I 

I 


Figure  Cl.  Sample  JH13189-4  13,  Delamination 
Surface  (200X) 


97 


Figure  C8 .  Sample  JH13189-4  #3,  Fracture 
Surface  (1500X) 
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Appendix  D:  Test  Results  from  LHMEL-1 


To  verify  the  plausibility  of  exploring  a  decrease  in 
effective  ablation  enthalpy  for  graphite/epoxy,  initial 
tests  were  conducted  using  the  15  kV  Laser  Hardened 
Materials  Evaluation  Laser-1  (LHMEL-1).  These  were  cursory 
tests  with  very  limited  objectives.  The  only  results 
intended  were  to  see  if  graphite/epoxy  did  Indeed  ablate 
faster  when  loaded  in  tension  than  when  unloaded.  These 
data  are  shown  in  Figure  Dl.  The  number  labeling  the  data 
points  are  intensities  at  which  each  was  taken.  Some  data 
on  temperature  rise  in  the  center  of  sample  ICH13989-4I12 
was  also  collected  (Figure  3).  This  sample  was  2.0  Inches 
(5.08  cm)  wide  and  0.244  inches  (0.620  cm)  thick  the 
thermocouple  plotted  in  Figure  3  was  in  the  center  of  the 
sample.  It  was  exposed  to  1700  kW/cm2  (11  kW  power  over 
6.53  cm2  area)  for  6.0  seconds.  These  parameters  were 
input  to  MELT 3D . FOR  to  generate  the  data  plotted  in  the 
figure . 


Table  III. 

LHMEL-1  Test  Results 


Power  Area 
(kW) (cm2) 

Intensity 

(kW/cm2) 

Time 

(sec) 

Mass 

Loss 

(gm) 

Q*ef  f 

(kJ/gm) 

LOAD 

( lbf )  (N) 

Stress 

Ratio 

11.5 

1.54 

7.5 

2.00 

0.65 

35.3 

12000 

53379 

0.33 

11.5 

1.54 

7.5 

1.19 

0.40 

34.2 

0 

0 

0.00 

11.5 

2.18 

5.3 

2.18 

0.75 

33.4 

12000 

53379 

0.50 

11.5 

2.18 

5.3 

3.00 

0.90 

38.3 

6000 

26689 

0.25 

11.5 

2.18 

5.3 

2.20 

0.60 

42.1 

9000 

40034 

0.37 

11.5 

2.18 

5.3 

6.50 

1.88 

39.7 

0 

0 

0.00 

99 


Q*  (kJ/gm) 


45 


o  5.3 

*0-1.5.* 


35 


30  _| - T - - — i i - 1 - 1 - 1 

0  0.1  0.2  0.3  0.4  0.5 

Stress  Ratio  (LOAD  /  LOAD  ujfimafe) 

Figure  D1 .  LHMEL-1  Preliminary  Results 
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